Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
lagrange_dirichlet_boundary_impl.h
Go to the documentation of this file.
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
39
40#ifdef UG_FOR_LUA
42#endif
43
44namespace ug{
45
46
48// setup
50
51template <typename TDomain, typename TAlgebra>
54{
55 base_type::set_approximation_space(approxSpace);
56 m_spApproxSpace = approxSpace;
57 m_spDomain = approxSpace->domain();
58 m_aaPos = m_spDomain->position_accessor();
59}
60
61template <typename TDomain, typename TAlgebra>
63clear()
64{
65 m_vBNDNumberData.clear();
66 m_vNumberData.clear();
67 m_vConstNumberData.clear();
68 m_vVectorData.clear();
69}
70
71template <typename TDomain, typename TAlgebra>
73add(SmartPtr<UserData<number, dim, bool> > func, const char* function, const char* subsets)
74{
75 m_vBNDNumberData.push_back(CondNumberData(func, function, subsets));
76}
77
78template <typename TDomain, typename TAlgebra>
80add(SmartPtr<UserData<number, dim, bool> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
81{
82 std::string function;
83 for(size_t i = 0; i < Fcts.size(); ++i){
84 if(i > 0) function.append(",");
85 function.append(Fcts[i]);
86 }
87 std::string subsets;
88 for(size_t i = 0; i < Subsets.size(); ++i){
89 if(i > 0) subsets.append(",");
90 subsets.append(Subsets[i]);
91 }
92
93 add(func, function.c_str(), subsets.c_str());
94}
95
96template <typename TDomain, typename TAlgebra>
98add(SmartPtr<UserData<number, dim> > func, const char* function, const char* subsets)
99{
100 m_vNumberData.push_back(NumberData(func, function, subsets));
101}
102
103template <typename TDomain, typename TAlgebra>
105add(SmartPtr<UserData<number, dim> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
106{
107 std::string function;
108 for(size_t i = 0; i < Fcts.size(); ++i){
109 if(i > 0) function.append(",");
110 function.append(Fcts[i]);
111 }
112 std::string subsets;
113 for(size_t i = 0; i < Subsets.size(); ++i){
114 if(i > 0) subsets.append(",");
115 subsets.append(Subsets[i]);
116 }
117
118 add(func, function.c_str(), subsets.c_str());
119}
120
121template <typename TDomain, typename TAlgebra>
123add(number value, const char* function, const char* subsets)
124{
125 m_vConstNumberData.push_back(ConstNumberData(value, function, subsets));
126}
127
128template <typename TDomain, typename TAlgebra>
130add(number value, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
131{
132 std::string function;
133 for(size_t i = 0; i < Fcts.size(); ++i){
134 if(i > 0) function.append(",");
135 function.append(Fcts[i]);
136 }
137 std::string subsets;
138 for(size_t i = 0; i < Subsets.size(); ++i){
139 if(i > 0) subsets.append(",");
140 subsets.append(Subsets[i]);
141 }
142
143 add(value, function.c_str(), subsets.c_str());
144}
145
146template <typename TDomain, typename TAlgebra>
148add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions, const char* subsets)
149{
150 m_vVectorData.push_back(VectorData(func, functions, subsets));
151}
152
153template <typename TDomain, typename TAlgebra>
155add(SmartPtr<UserData<MathVector<dim>, dim> > func, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
156{
157 std::string function;
158 for(size_t i = 0; i < Fcts.size(); ++i){
159 if(i > 0) function.append(",");
160 function.append(Fcts[i]);
161 }
162 std::string subsets;
163 for(size_t i = 0; i < Subsets.size(); ++i){
164 if(i > 0) subsets.append(",");
165 subsets.append(Subsets[i]);
166 }
167
168 add(func, function.c_str(), subsets.c_str());
169}
170
171#ifdef UG_FOR_LUA
172template <typename TDomain, typename TAlgebra>
174add(const char* name, const char* function, const char* subsets)
175{
179 add(sp, function, subsets);
180 return;
181 }
185 add(sp, function, subsets);
186 return;
187 }
188 if(LuaUserData<MathVector<dim>, dim>::check_callback_returns(name)){
190 LuaUserDataFactory<MathVector<dim>, dim>::create(name);
191 add(sp, function, subsets);
192 return;
193 }
194
195// no match found
197 UG_THROW("LagrangeDirichlet::add: Lua-Callback with name '"<<name<<
198 "' does not exist.");
199
200// name exists but wrong signature
201 UG_THROW("LagrangeDirichlet::add: Cannot find matching callback "
202 "signature. Use one of:\n"
203 "a) Number - Callback\n"
205 "b) Conditional Number - Callback\n"
207 "c) "<<dim<<"d Vector - Callback\n"
208 << (LuaUserData<MathVector<dim>, dim>::signature()));
209}
210
211template <typename TDomain, typename TAlgebra>
213add(const char* name, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
214{
215 std::string function;
216 for(size_t i = 0; i < Fcts.size(); ++i){
217 if(i > 0) function.append(",");
218 function.append(Fcts[i]);
219 }
220 std::string subsets;
221 for(size_t i = 0; i < Subsets.size(); ++i){
222 if(i > 0) subsets.append(",");
223 subsets.append(Subsets[i]);
224 }
225
226 add(name, function.c_str(), subsets.c_str());
227}
228
229template <typename TDomain, typename TAlgebra>
231add(LuaFunctionHandle fct, const char* function, const char* subsets)
232{
235 make_sp(new LuaUserData<number, dim>(fct));
236 add(sp, function, subsets);
237 return;
238 }
241 make_sp(new LuaUserData<number, dim, bool>(fct));
242 add(sp, function, subsets);
243 return;
244 }
245 if(LuaUserData<MathVector<dim>, dim>::check_callback_returns(fct)){
247 make_sp(new LuaUserData<MathVector<dim>, dim>(fct));
248 add(sp, function, subsets);
249 return;
250 }
251
252// name exists but wrong signature
253 UG_THROW("LagrangeDirichlet::add: Cannot find matching callback "
254 "signature. Use one of:\n"
255 "a) Number - Callback\n"
257 "b) Conditional Number - Callback\n"
259 "c) "<<dim<<"d Vector - Callback\n"
260 << (LuaUserData<MathVector<dim>, dim>::signature()));
261}
262
263template <typename TDomain, typename TAlgebra>
265add(LuaFunctionHandle fct, const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
266{
267 std::string function;
268 for(size_t i = 0; i < Fcts.size(); ++i){
269 if(i > 0) function.append(",");
270 function.append(Fcts[i]);
271 }
272 std::string subsets;
273 for(size_t i = 0; i < Subsets.size(); ++i){
274 if(i > 0) subsets.append(",");
275 subsets.append(Subsets[i]);
276 }
277
278 add(fct, function.c_str(), subsets.c_str());
279}
280#endif
281
282template <typename TDomain, typename TAlgebra>
284add(const char* functions, const char* subsets)
285{
286 m_vOldNumberData.push_back(OldNumberData(functions, subsets));
287}
288
289template <typename TDomain, typename TAlgebra>
291add(const std::vector<std::string>& Fcts, const std::vector<std::string>& Subsets)
292{
293 std::string function;
294 for(size_t i = 0; i < Fcts.size(); ++i){
295 if(i > 0) function.append(",");
296 function.append(Fcts[i]);
297 }
298 std::string subsets;
299 for(size_t i = 0; i < Subsets.size(); ++i){
300 if(i > 0) subsets.append(",");
301 subsets.append(Subsets[i]);
302 }
303
304 add(function.c_str(), subsets.c_str());
305}
306
307template <typename TDomain, typename TAlgebra>
309check_functions_and_subsets(FunctionGroup& functionGroup, SubsetGroup& subsetGroup, size_t numFct) const
310{
311// only number of functions allowed
312 if(functionGroup.size() != numFct)
313 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 ConstSmartPtr<ISubsetHandler> pSH = m_spApproxSpace->subset_handler();
320
321// loop subsets
322 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 if(subsetIndex < 0 || subsetIndex >= pSH->num_subsets())
329 UG_THROW("DirichletBoundary:extract_data:"
330 " Invalid Subset Index " << subsetIndex << ". (Valid is"
331 " 0, .. , " << pSH->num_subsets() <<").");
332
333 // check all functions
334 for(size_t i=0; i < functionGroup.size(); ++i)
335 {
336 const size_t fct = functionGroup[i];
337
338 // check if function exist
339 if(fct >= m_spApproxSpace->num_fct())
340 UG_THROW("DirichletBoundary:extract_data:"
341 " Function "<< fct << " does not exist in pattern.");
342
343 // check that function is defined for segment
344 if(!m_spApproxSpace->is_def_in_subset(fct, subsetIndex))
345 UG_THROW("DirichletBoundary:extract_data:"
346 " Function "<<fct<<" not defined on subset "<<subsetIndex);
347 }
348 }
349
350// everything ok
351}
352
353template <typename TDomain, typename TAlgebra>
354template <typename TUserData, typename TScheduledUserData>
356extract_data(std::map<int, std::vector<TUserData*> >& mvUserDataBndSegment,
357 std::vector<TScheduledUserData>& vUserData)
358{
359// clear the extracted data
360 mvUserDataBndSegment.clear();
361
362 for(size_t i = 0; i < vUserData.size(); ++i)
363 {
364 // create Function Group and Subset Group
365 try
366 {
367 if (! m_bInvertSubsetSelection)
368 vUserData[i].ssGrp = m_spApproxSpace->subset_grp_by_name
369 (vUserData[i].ssName.c_str());
370 else
371 vUserData[i].ssGrp = m_spApproxSpace->all_subsets_grp_except_for
372 (vUserData[i].ssName.c_str());
373 }
374 UG_CATCH_THROW(" Subsets '"<<vUserData[i].ssName<<"' not"
375 " all contained in ApproximationSpace.");
376
377 FunctionGroup fctGrp;
378 try{
379 fctGrp = m_spApproxSpace->fct_grp_by_name(vUserData[i].fctName.c_str());
380 }UG_CATCH_THROW(" Functions '"<<vUserData[i].fctName<<"' not"
381 " all contained in ApproximationSpace.");
382
383 // check functions and subsets
384 check_functions_and_subsets(fctGrp, vUserData[i].ssGrp, TUserData::numFct);
385
386 // set functions
387 if(fctGrp.size() != TUserData::numFct)
388 UG_THROW("LagrangeDirichletBoundary: wrong number of fct");
389
390 for(size_t fct = 0; fct < TUserData::numFct; ++fct)
391 {
392 vUserData[i].fct[fct] = fctGrp[fct];
393 }
394
395 // loop subsets
396 for(size_t si = 0; si < vUserData[i].ssGrp.size(); ++si)
397 {
398 // get subset index and function
399 const int subsetIndex = vUserData[i].ssGrp[si];
400
401 // remember functor and function
402 mvUserDataBndSegment[subsetIndex].push_back(&vUserData[i]);
403 }
404 }
405}
406
407
408template <typename TDomain, typename TAlgebra>
411{
412// check that function pattern exists
413 if(!m_spApproxSpace.valid())
414 UG_THROW("DirichletBoundary:extract_data: "
415 " Approximation Space not set.");
416
417 extract_data(m_mNumberBndSegment, m_vNumberData);
418 extract_data(m_mBNDNumberBndSegment, m_vBNDNumberData);
419 extract_data(m_mConstNumberBndSegment, m_vConstNumberData);
420 extract_data(m_mVectorBndSegment, m_vVectorData);
421 extract_data(m_mOldNumberBndSegment, m_vOldNumberData);
422}
423
425// assemble_dirichlet_rows
427
428template <typename TDomain, typename TAlgebra>
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
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
480// adjust TRANSFER
482
483
484template <typename TDomain, typename TAlgebra>
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 extract_data();
500
501 adjust_prolongation<CondNumberData>(m_mBNDNumberBndSegment, P, ddFine, ddCoarse, time);
502 adjust_prolongation<NumberData>(m_mNumberBndSegment, P, ddFine, ddCoarse, time);
503 adjust_prolongation<ConstNumberData>(m_mConstNumberBndSegment, P, ddFine, ddCoarse, time);
504
505 adjust_prolongation<VectorData>(m_mVectorBndSegment, P, ddFine, ddCoarse, time);
506
507 adjust_prolongation<OldNumberData>(m_mOldNumberBndSegment, P, ddFine, ddCoarse, time);
508}
509
510template <typename TDomain, typename TAlgebra>
511template <typename TUserData>
513adjust_prolongation(const std::map<int, std::vector<TUserData*> >& mvUserData,
514 matrix_type& P,
517 number time)
518{
519// loop boundary subsets
520 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
521 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
522 {
523 // get subset index
524 const int si = (*iter).first;
525
526 // get vector of scheduled dirichlet data on this subset
527 const std::vector<TUserData*>& vUserData = (*iter).second;
528
529 // adapt jacobian for dofs in each base element type
530 try
531 {
532 if(ddFine->max_dofs(VERTEX)) adjust_prolongation<RegularVertex, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
533 if(ddFine->max_dofs(EDGE)) adjust_prolongation<Edge, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
534 if(ddFine->max_dofs(FACE)) adjust_prolongation<Face, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
535 if(ddFine->max_dofs(VOLUME)) adjust_prolongation<Volume, TUserData>(vUserData, si, P, ddFine, ddCoarse, time);
536 }
537 UG_CATCH_THROW("DirichletBoundary::adjust_prolongation:"
538 " While calling 'adapt_prolongation' for TUserData, aborting.");
539 }
540}
541
542template <typename TDomain, typename TAlgebra>
543template <typename TBaseElem, typename TUserData>
545adjust_prolongation(const std::vector<TUserData*>& vUserData, int si,
546 matrix_type& P,
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
562 iter = ddFine->begin<TBaseElem>(si);
563 iterEnd = ddFine->end<TBaseElem>(si);
564
565// loop elements
566 for( ; iter != iterEnd; iter++)
567 {
568 // get vertex
569 TBaseElem* elem = *iter;
570 GridObject* parent = m_spDomain->grid()->get_parent(elem);
571 if(!parent) continue;
572 if(!ddCoarse->is_contained(parent)) continue;
573
574 // loop dirichlet functions on this segment
575 for(size_t i = 0; i < vUserData.size(); ++i)
576 {
577 for(size_t f = 0; f < TUserData::numFct; ++f)
578 {
579 // get function index
580 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 ddFine->inner_dof_indices(elem, fct, vFineDoF);
587 ddCoarse->inner_dof_indices(parent, fct, vCoarseDoF);
588
589 // get dof position
590 if(TUserData::isConditional){
591 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
592 UG_ASSERT(vFineDoF.size() == vPos.size(), "Size mismatch");
593 }
594
595 // loop dofs on element
596 for(size_t j = 0; j < vFineDoF.size(); ++j)
597 {
598 // check if function is dirichlet
599 if(TUserData::isConditional){
600 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
601 }
602
603 SetRow(P, vFineDoF[j], 0.0);
604 }
605
606 if(vFineDoF.size() > 0){
607 for(size_t k = 0; k < vCoarseDoF.size(); ++k){
608 DoFRef(P, vFineDoF[0], vCoarseDoF[k]) = 1.0;
609 }
610 }
611 }
612 }
613 }
614}
615
616template <typename TDomain, typename TAlgebra>
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 extract_data();
632
633 adjust_restriction<CondNumberData>(m_mBNDNumberBndSegment, R, ddCoarse, ddFine, time);
634 adjust_restriction<NumberData>(m_mNumberBndSegment, R, ddCoarse, ddFine, time);
635 adjust_restriction<ConstNumberData>(m_mConstNumberBndSegment, R, ddCoarse, ddFine, time);
636
637 adjust_restriction<VectorData>(m_mVectorBndSegment, R, ddCoarse, ddFine, time);
638
639 adjust_restriction<OldNumberData>(m_mOldNumberBndSegment, R, ddCoarse, ddFine, time);
640}
641
642template <typename TDomain, typename TAlgebra>
643template <typename TUserData>
645adjust_restriction(const std::map<int, std::vector<TUserData*> >& mvUserData,
646 matrix_type& R,
649 number time)
650{
651// loop boundary subsets
652 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
653 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
654 {
655 // get subset index
656 const int si = (*iter).first;
657
658 // get vector of scheduled dirichlet data on this subset
659 const std::vector<TUserData*>& vUserData = (*iter).second;
660
661 // adapt jacobian for dofs in each base element type
662 try
663 {
664 if(ddFine->max_dofs(VERTEX)) adjust_restriction<RegularVertex, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
665 if(ddFine->max_dofs(EDGE)) adjust_restriction<Edge, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
666 if(ddFine->max_dofs(FACE)) adjust_restriction<Face, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
667 if(ddFine->max_dofs(VOLUME)) adjust_restriction<Volume, TUserData>(vUserData, si, R, ddCoarse, ddFine, time);
668 }
669 UG_CATCH_THROW("DirichletBoundary::adjust_restriction:"
670 " While calling 'adjust_restriction' for TUserData, aborting.");
671 }
672}
673
674template <typename TDomain, typename TAlgebra>
675template <typename TBaseElem, typename TUserData>
677adjust_restriction(const std::vector<TUserData*>& vUserData, int si,
678 matrix_type& R,
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
694 iter = ddFine->begin<TBaseElem>(si);
695 iterEnd = ddFine->end<TBaseElem>(si);
696
697// loop elements
698 for( ; iter != iterEnd; iter++)
699 {
700 // get vertex
701 TBaseElem* elem = *iter;
702 GridObject* parent = m_spDomain->grid()->get_parent(elem);
703 if(!parent) continue;
704 if(!ddCoarse->is_contained(parent)) continue;
705
706 // loop dirichlet functions on this segment
707 for(size_t i = 0; i < vUserData.size(); ++i)
708 {
709 for(size_t f = 0; f < TUserData::numFct; ++f)
710 {
711 // get function index
712 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 ddFine->inner_dof_indices(elem, fct, vFineDoF);
719 ddCoarse->inner_dof_indices(parent, fct, vCoarseDoF);
720
721 // get dof position
722 if(TUserData::isConditional){
723 InnerDoFPosition<TDomain>(vPos, parent, *m_spDomain, lfeID);
724 UG_ASSERT(vCoarseDoF.size() == vPos.size(), "Size mismatch");
725 }
726
727 // loop dofs on element
728 for(size_t j = 0; j < vCoarseDoF.size(); ++j)
729 {
730 // check if function is dirichlet
731 if(TUserData::isConditional){
732 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
733 }
734
735 SetRow(R, vCoarseDoF[j], 0.0);
736 }
737
738 if(vFineDoF.size() > 0){
739 for(size_t k = 0; k < vCoarseDoF.size(); ++k){
740 DoFRef(R, vCoarseDoF[k], vFineDoF[0]) = 1.0;
741 }
742 }
743 }
744 }
745 }
746}
747
749// adjust JACOBIAN
751
752template <typename TDomain, typename TAlgebra>
755 ConstSmartPtr<DoFDistribution> dd, int type, number time,
757 const number s_a0)
758{
759 extract_data();
760
761 adjust_jacobian<CondNumberData>(m_mBNDNumberBndSegment, J, u, dd, time);
762 adjust_jacobian<NumberData>(m_mNumberBndSegment, J, u, dd, time);
763 adjust_jacobian<ConstNumberData>(m_mConstNumberBndSegment, J, u, dd, time);
764
765 adjust_jacobian<VectorData>(m_mVectorBndSegment, J, u, dd, time);
766
767 adjust_jacobian<OldNumberData>(m_mOldNumberBndSegment, J, u, dd, time);
768}
769
770
771template <typename TDomain, typename TAlgebra>
772template <typename TUserData>
774adjust_jacobian(const std::map<int, std::vector<TUserData*> >& mvUserData,
775 matrix_type& J, const vector_type& u,
777{
778// loop boundary subsets
779 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
780 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
781 {
782 // get subset index
783 const int si = (*iter).first;
784
785 // get vector of scheduled dirichlet data on this subset
786 const std::vector<TUserData*>& vUserData = (*iter).second;
787
788 // adapt jacobian for dofs in each base element type
789 try
790 {
791 if(dd->max_dofs(VERTEX))
792 adjust_jacobian<RegularVertex, TUserData>(vUserData, si, J, u, dd, time);
793 if(dd->max_dofs(EDGE))
794 adjust_jacobian<Edge, TUserData>(vUserData, si, J, u, dd, time);
795 if(dd->max_dofs(FACE))
796 adjust_jacobian<Face, TUserData>(vUserData, si, J, u, dd, time);
797 if(dd->max_dofs(VOLUME))
798 adjust_jacobian<Volume, TUserData>(vUserData, si, J, u, dd, time);
799 }
800 UG_CATCH_THROW("DirichletBoundary::adjust_jacobian:"
801 " While calling 'adapt_jacobian' for TUserData, aborting.");
802 }
803}
804
805template <typename TDomain, typename TAlgebra>
806template <typename TBaseElem, typename TUserData>
808adjust_jacobian(const std::vector<TUserData*>& vUserData, int si,
809 matrix_type& J, const vector_type& u,
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
827 iter = dd->begin<TBaseElem>(si);
828 iterEnd = dd->end<TBaseElem>(si);
829
830// loop elements
831 for( ; iter != iterEnd; iter++)
832 {
833 // get vertex
834 TBaseElem* elem = *iter;
835
836 // loop dirichlet functions on this segment
837 for(size_t i = 0; i < vUserData.size(); ++i)
838 {
839 for(size_t f = 0; f < TUserData::numFct; ++f)
840 {
841 // get function index
842 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 dd->inner_dof_indices(elem, fct, multInd);
849
850 // get dof position
851 if(TUserData::isConditional){
852 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
853 UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
854 }
855
856 // loop dofs on element
857 for(size_t j = 0; j < multInd.size(); ++j)
858 {
859 // check if function is dirichlet
860 if(TUserData::isConditional){
861 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
862 }
863
864 this->m_spAssTuner->set_dirichlet_row(J, multInd[j]);
865 if(m_bDirichletColumns)
866 dirichletDoFIndices.insert(multInd[j][0]);
867 }
868 }
869 }
870 }
871
872
873 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 for(size_t i = 0; i<nr; i++)
885 {
886 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 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 if(i!=it.index())
905 it.value() = 0.0;
906 }
907
908 }
909 }
910 }
911
912}
913
915// adjust DEFECT
917
918template <typename TDomain, typename TAlgebra>
921 ConstSmartPtr<DoFDistribution> dd, int type, number time,
923 const std::vector<number>* vScaleMass,
924 const std::vector<number>* vScaleStiff)
925{
926 extract_data();
927
928 adjust_defect<CondNumberData>(m_mBNDNumberBndSegment, d, u, dd, time);
929 adjust_defect<NumberData>(m_mNumberBndSegment, d, u, dd, time);
930 adjust_defect<ConstNumberData>(m_mConstNumberBndSegment, d, u, dd, time);
931
932 adjust_defect<VectorData>(m_mVectorBndSegment, d, u, dd, time);
933
934 adjust_defect<OldNumberData>(m_mOldNumberBndSegment, d, u, dd, time);
935}
936
937
938template <typename TDomain, typename TAlgebra>
939template <typename TUserData>
941adjust_defect(const std::map<int, std::vector<TUserData*> >& mvUserData,
942 vector_type& d, const vector_type& u,
944{
945// loop boundary subsets
946 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
947 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
948 {
949 // get subset index
950 const int si = (*iter).first;
951
952 // get vector of scheduled dirichlet data on this subset
953 const std::vector<TUserData*>& vUserData = (*iter).second;
954
955 // adapt jacobian for dofs in each base element type
956 try
957 {
958 if(dd->max_dofs(VERTEX))
959 adjust_defect<RegularVertex, TUserData>(vUserData, si, d, u, dd, time);
960 if(dd->max_dofs(EDGE))
961 adjust_defect<Edge, TUserData>(vUserData, si, d, u, dd, time);
962 if(dd->max_dofs(FACE))
963 adjust_defect<Face, TUserData>(vUserData, si, d, u, dd, time);
964 if(dd->max_dofs(VOLUME))
965 adjust_defect<Volume, TUserData>(vUserData, si, d, u, dd, time);
966 }
967 UG_CATCH_THROW("DirichletBoundary::adjust_defect:"
968 " While calling 'adjust_defect' for TUserData, aborting.");
969 }
970}
971
972template <typename TDomain, typename TAlgebra>
973template <typename TBaseElem, typename TUserData>
975adjust_defect(const std::vector<TUserData*>& vUserData, int si,
976 vector_type& d, const vector_type& u,
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
990 iter = dd->begin<TBaseElem>(si);
991 iterEnd = dd->end<TBaseElem>(si);
992
993// loop elements
994 for( ; iter != iterEnd; iter++)
995 {
996 // get vertex
997 TBaseElem* elem = *iter;
998
999 // loop dirichlet functions on this segment
1000 for(size_t i = 0; i < vUserData.size(); ++i)
1001 {
1002 for(size_t f = 0; f < TUserData::numFct; ++f)
1003 {
1004 // get function index
1005 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 dd->inner_dof_indices(elem, fct, multInd);
1012
1013 // get dof position
1014 if(TUserData::isConditional){
1015 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 for(size_t j = 0; j < multInd.size(); ++j)
1022 {
1023 // check if function is dirichlet
1024 if(TUserData::isConditional){
1025 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1026 }
1027
1028 // set zero for dirichlet values
1029 this->m_spAssTuner->set_dirichlet_val(d, multInd[j], 0.0);
1030 }
1031 }
1032 }
1033 }
1034}
1035
1037// adjust SOLUTION
1039
1040template <typename TDomain, typename TAlgebra>
1043{
1044 extract_data();
1045
1046 adjust_solution<CondNumberData>(m_mBNDNumberBndSegment, u, dd, time);
1047 adjust_solution<NumberData>(m_mNumberBndSegment, u, dd, time);
1048 adjust_solution<ConstNumberData>(m_mConstNumberBndSegment, u, dd, time);
1049
1050 adjust_solution<VectorData>(m_mVectorBndSegment, u, dd, time);
1051
1052 adjust_solution<OldNumberData>(m_mOldNumberBndSegment, u, dd, time);
1053}
1054
1055
1056template <typename TDomain, typename TAlgebra>
1057template <typename TUserData>
1059adjust_solution(const std::map<int, std::vector<TUserData*> >& mvUserData,
1061{
1062// loop boundary subsets
1063 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1064 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1065 {
1066 // get subset index
1067 const int si = (*iter).first;
1068
1069 // get vector of scheduled dirichlet data on this subset
1070 const std::vector<TUserData*>& vUserData = (*iter).second;
1071
1072 // adapt jacobian for dofs in each base element type
1073 try
1074 {
1075 if(dd->max_dofs(VERTEX))
1076 adjust_solution<RegularVertex, TUserData>(vUserData, si, u, dd, time);
1077 if(dd->max_dofs(EDGE))
1078 adjust_solution<Edge, TUserData>(vUserData, si, u, dd, time);
1079 if(dd->max_dofs(FACE))
1080 adjust_solution<Face, TUserData>(vUserData, si, u, dd, time);
1081 if(dd->max_dofs(VOLUME))
1082 adjust_solution<Volume, TUserData>(vUserData, si, u, dd, time);
1083 }
1084 UG_CATCH_THROW("DirichletBoundary::adjust_solution:"
1085 " While calling 'adjust_solution' for TUserData, aborting.");
1086 }
1087}
1088
1089template <typename TDomain, typename TAlgebra>
1090template <typename TBaseElem, typename TUserData>
1092adjust_solution(const std::vector<TUserData*>& vUserData, int si,
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
1110 iter = dd->begin<TBaseElem>(si);
1111 iterEnd = dd->end<TBaseElem>(si);
1112
1113// loop elements
1114 for( ; iter != iterEnd; iter++)
1115 {
1116 // get vertex
1117 TBaseElem* elem = *iter;
1118
1119 // loop dirichlet functions on this segment
1120 for(size_t i = 0; i < vUserData.size(); ++i)
1121 {
1122 for(size_t f = 0; f < TUserData::numFct; ++f)
1123 {
1124 // get function index
1125 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 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1132
1133 // get multi indices
1134 dd->inner_dof_indices(elem, fct, multInd);
1135
1136 UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
1137
1138 // loop dofs on element
1139 for(size_t j = 0; j < multInd.size(); ++j)
1140 {
1141 // get dirichlet value
1142 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1143
1144 this->m_spAssTuner->set_dirichlet_val(u, multInd[j], val[f]);
1145 }
1146 }
1147 }
1148 }
1149}
1150
1151
1152
1154// adjust CORRECTION
1156
1157template <typename TDomain, typename TAlgebra>
1159(
1160 vector_type& c,
1162 int type,
1163 number time
1164)
1165{
1166 extract_data();
1167
1168 adjust_correction<CondNumberData>(m_mBNDNumberBndSegment, c, dd, time);
1169 adjust_correction<NumberData>(m_mNumberBndSegment, c, dd, time);
1170 adjust_correction<ConstNumberData>(m_mConstNumberBndSegment, c, dd, time);
1171
1172 adjust_correction<VectorData>(m_mVectorBndSegment, c, dd, time);
1173
1174 adjust_correction<OldNumberData>(m_mOldNumberBndSegment, c, dd, time);
1175}
1176
1177template <typename TDomain, typename TAlgebra>
1178template <typename TUserData>
1180adjust_correction(const std::map<int, std::vector<TUserData*> >& mvUserData,
1182{
1183// loop boundary subsets
1184 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1185 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1186 {
1187 // get subset index
1188 const int si = (*iter).first;
1189
1190 // get vector of scheduled dirichlet data on this subset
1191 const std::vector<TUserData*>& vUserData = (*iter).second;
1192
1193 // adapt correction for dofs in each base element type
1194 try
1195 {
1196 if(dd->max_dofs(VERTEX))
1197 adjust_correction<RegularVertex, TUserData>(vUserData, si, c, dd, time);
1198 if(dd->max_dofs(EDGE))
1199 adjust_correction<Edge, TUserData>(vUserData, si,c, dd, time);
1200 if(dd->max_dofs(FACE))
1201 adjust_correction<Face, TUserData>(vUserData, si, c, dd, time);
1202 if(dd->max_dofs(VOLUME))
1203 adjust_correction<Volume, TUserData>(vUserData, si, c, dd, time);
1204 }
1205 UG_CATCH_THROW("DirichletBoundary::adjust_correction:"
1206 " While calling 'adjust_correction' for TUserData, aborting.");
1207 }
1208}
1209
1210template <typename TDomain, typename TAlgebra>
1211template <typename TBaseElem, typename TUserData>
1213adjust_correction(const std::vector<TUserData*>& vUserData, int si,
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
1227 iter = dd->begin<TBaseElem>(si);
1228 iterEnd = dd->end<TBaseElem>(si);
1229
1230// loop elements
1231 for( ; iter != iterEnd; iter++)
1232 {
1233 // get vertex
1234 TBaseElem* elem = *iter;
1235
1236 // loop dirichlet functions on this segment
1237 for(size_t i = 0; i < vUserData.size(); ++i)
1238 {
1239 for(size_t f = 0; f < TUserData::numFct; ++f)
1240 {
1241 // get function index
1242 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 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1249
1250 // get multi indices
1251 dd->inner_dof_indices(elem, fct, multInd);
1252
1253 UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
1254
1255 // loop dofs on element
1256 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 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1260
1261 this->m_spAssTuner->set_dirichlet_val(c, multInd[j], 0.0);
1262 }
1263 }
1264 }
1265 }
1266}
1267
1268
1270// adjust LINEAR
1272
1273template <typename TDomain, typename TAlgebra>
1276 ConstSmartPtr<DoFDistribution> dd, int type, number time)
1277{
1278 extract_data();
1279
1280 adjust_linear<CondNumberData>(m_mBNDNumberBndSegment, A, b, dd, time);
1281 adjust_linear<NumberData>(m_mNumberBndSegment, A, b, dd, time);
1282 adjust_linear<ConstNumberData>(m_mConstNumberBndSegment, A, b, dd, time);
1283
1284 adjust_linear<VectorData>(m_mVectorBndSegment, A, b, dd, time);
1285
1286 adjust_linear<OldNumberData>(m_mOldNumberBndSegment, A, b, dd, time);
1287}
1288
1289
1290template <typename TDomain, typename TAlgebra>
1291template <typename TUserData>
1293adjust_linear(const std::map<int, std::vector<TUserData*> >& mvUserData,
1294 matrix_type& A, vector_type& b,
1296{
1297// loop boundary subsets
1298 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1299 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1300 {
1301 // get subset index
1302 const int si = (*iter).first;
1303
1304 // get vector of scheduled dirichlet data on this subset
1305 const std::vector<TUserData*>& vUserData = (*iter).second;
1306
1307 // adapt jacobian for dofs in each base element type
1308 try
1309 {
1310 if(dd->max_dofs(VERTEX))
1311 adjust_linear<RegularVertex, TUserData>(vUserData, si, A, b, dd, time);
1312 if(dd->max_dofs(EDGE))
1313 adjust_linear<Edge, TUserData>(vUserData, si, A, b, dd, time);
1314 if(dd->max_dofs(FACE))
1315 adjust_linear<Face, TUserData>(vUserData, si, A, b, dd, time);
1316 if(dd->max_dofs(VOLUME))
1317 adjust_linear<Volume, TUserData>(vUserData, si, A, b, dd, time);
1318 }
1319 UG_CATCH_THROW("DirichletBoundary::adjust_linear:"
1320 " While calling 'adjust_linear' for TUserData, aborting.");
1321 }
1322}
1323
1324template <typename TDomain, typename TAlgebra>
1325template <typename TBaseElem, typename TUserData>
1327adjust_linear(const std::vector<TUserData*>& vUserData, int si,
1328 matrix_type& A, vector_type& b,
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
1345 iter = dd->begin<TBaseElem>(si);
1346 iterEnd = dd->end<TBaseElem>(si);
1347
1348// loop elements
1349 for( ; iter != iterEnd; iter++)
1350 {
1351 // get vertex
1352 TBaseElem* elem = *iter;
1353
1354 // loop dirichlet functions on this segment
1355 for(size_t i = 0; i < vUserData.size(); ++i)
1356 {
1357 for(size_t f = 0; f < TUserData::numFct; ++f)
1358 {
1359 // get function index
1360 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 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1367
1368 // get multi indices
1369 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 for(size_t j = 0; j < multInd.size(); ++j)
1377 {
1378 // check if function is dirichlet and read value
1379 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1380
1381 this->m_spAssTuner->set_dirichlet_row(A, multInd[j]);
1382
1383 if(m_bDirichletColumns)
1384 {
1385 // FIXME: Beware, this is dangerous!
1386 // It will not work for blocked algebras.
1387 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 this->m_spAssTuner->set_dirichlet_val(b, multInd[j], val[f]);
1394 }
1395 }
1396 }
1397 }
1398
1399
1400
1401 if(m_bDirichletColumns){
1402// UG_LOG("adjust linear\n")
1403 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 for(size_t i = 0; i<nr; i++)
1413 {
1414 // do not save any entries in Dirichlet rows!
1415 if (dirichletDoFIndices.find(i) != dirichletDoFIndices.end())
1416 continue;
1417
1418 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 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 m_dirichletMap[si][i][it.index()] = it.value();
1428 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 for(size_t i = 0; i<nr; i++)
1437 {
1438 // step over if this row did not contain any connections to Dirichlet nodes
1439 if ((itdirichletMap = m_dirichletMap[si].find(i)) == m_dirichletMap[si].end())
1440 continue;
1441
1442 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 if ((itdirichletRowMap = itdirichletMap->second.find(it.index())) != itdirichletMap->second.end())
1446 b[i] -= itdirichletRowMap->second*b[it.index()];
1447 }
1448 }
1449 }
1450}
1451
1453// adjust RHS
1455
1456template <typename TDomain, typename TAlgebra>
1459 ConstSmartPtr<DoFDistribution> dd, int type, number time)
1460{
1461 extract_data();
1462
1463 adjust_rhs<CondNumberData>(m_mBNDNumberBndSegment, b, u, dd, time);
1464 adjust_rhs<NumberData>(m_mNumberBndSegment, b, u, dd, time);
1465 adjust_rhs<ConstNumberData>(m_mConstNumberBndSegment, b, u, dd, time);
1466
1467 adjust_rhs<VectorData>(m_mVectorBndSegment, b, u, dd, time);
1468
1469 adjust_rhs<OldNumberData>(m_mOldNumberBndSegment, b, u, dd, time);
1470}
1471
1472
1473template <typename TDomain, typename TAlgebra>
1474template <typename TUserData>
1476adjust_rhs(const std::map<int, std::vector<TUserData*> >& mvUserData,
1477 vector_type& b, const vector_type& u,
1479{
1480// loop boundary subsets
1481 typename std::map<int, std::vector<TUserData*> >::const_iterator iter;
1482 for(iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1483 {
1484 // get subset index
1485 const int si = (*iter).first;
1486
1487 // get vector of scheduled dirichlet data on this subset
1488 const std::vector<TUserData*>& vUserData = (*iter).second;
1489
1490 // adapt jacobian for dofs in each base element type
1491 try
1492 {
1493 if(dd->max_dofs(VERTEX))
1494 adjust_rhs<RegularVertex, TUserData>(vUserData, si, b, u, dd, time);
1495 if(dd->max_dofs(EDGE))
1496 adjust_rhs<Edge, TUserData>(vUserData, si, b, u, dd, time);
1497 if(dd->max_dofs(FACE))
1498 adjust_rhs<Face, TUserData>(vUserData, si, b, u, dd, time);
1499 if(dd->max_dofs(VOLUME))
1500 adjust_rhs<Volume, TUserData>(vUserData, si, b, u, dd, time);
1501 }
1502 UG_CATCH_THROW("DirichletBoundary::adjust_rhs:"
1503 " While calling 'adjust_rhs' for TUserData, aborting.");
1504 }
1505}
1506
1507template <typename TDomain, typename TAlgebra>
1508template <typename TBaseElem, typename TUserData>
1510adjust_rhs(const std::vector<TUserData*>& vUserData, int si,
1511 vector_type& b, const vector_type& u,
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
1525 iter = dd->begin<TBaseElem>(si);
1526 iterEnd = dd->end<TBaseElem>(si);
1527
1528// loop elements
1529 for( ; iter != iterEnd; iter++)
1530 {
1531 // get vertex
1532 TBaseElem* elem = *iter;
1533
1534 // loop dirichlet functions on this segment
1535 for(size_t i = 0; i < vUserData.size(); ++i)
1536 {
1537 for(size_t f = 0; f < TUserData::numFct; ++f)
1538 {
1539 // get function index
1540 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 InnerDoFPosition<TDomain>(vPos, elem, *m_spDomain, lfeID);
1547
1548 // get multi indices
1549 dd->inner_dof_indices(elem, fct, multInd);
1550
1551 UG_ASSERT(multInd.size() == vPos.size(), "Size mismatch");
1552
1553 // loop dofs on element
1554 for(size_t j = 0; j < multInd.size(); ++j)
1555 {
1556 // check if function is dirichlet and read value
1557 if(!(*vUserData[i])(val, vPos[j], time, si)) continue;
1558
1559 if (TUserData::setSolValue)
1560 this->m_spAssTuner->set_dirichlet_val(b, multInd[j], val[f]);
1561 else
1562 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 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 size_t nr = m_A->num_rows();
1575 for(size_t i = 0; i<nr; i++)
1576 {
1577 // step over if this row did not contain any connections to Dirichlet nodes
1578 if ((itdirichletMap = m_dirichletMap[si].find(i)) == m_dirichletMap[si].end())
1579 continue;
1580
1581 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 if ((itdirichletRowMap = itdirichletMap->second.find(it.index())) != itdirichletMap->second.end())
1585 b[i] -= itdirichletRowMap->second*b[it.index()];
1586 }
1587 }
1588 }
1589}
1590
1591
1592// //////////////////////////////////////////////////////////////////////////////
1593// adjust error
1594// //////////////////////////////////////////////////////////////////////////////
1595
1596template <typename TDomain, typename TAlgebra>
1599(
1600 const vector_type& u,
1602 int type,
1603 number time,
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 if (this->m_spErrEstData.get() == NULL)
1611 {
1612 UG_THROW("No ErrEstData object has been given to this constraint!");
1613 }
1614
1615 err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
1616
1617 if (!err_est_data)
1618 {
1619 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 extract_data();
1625
1626 adjust_error<CondNumberData>(m_mBNDNumberBndSegment, u, dd, time);
1627 adjust_error<NumberData>(m_mNumberBndSegment, u, dd, time);
1628 adjust_error<ConstNumberData>(m_mConstNumberBndSegment, u, dd, time);
1629
1630 adjust_error<VectorData>(m_mVectorBndSegment, u, dd, time);
1631
1632 adjust_error<OldNumberData>(m_mOldNumberBndSegment, u, dd, time);
1633}
1634
1635template <typename TDomain, typename TAlgebra>
1636template <typename TUserData>
1639(
1640 const std::map<int, std::vector<TUserData*> >& mvUserData,
1641 const vector_type& u,
1643 number time
1644)
1645{
1646 // cast error estimator data object to the right type
1647 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 for (iter = mvUserData.begin(); iter != mvUserData.end(); ++iter)
1654 {
1655 // get subset index
1656 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
1674 iter = dd->begin<side_type>(si);
1675 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 for( ; iter != iterEnd; iter++)
1680 {
1681 // get vertex
1682 side_type* elem = *iter;
1683
1684 // get reference object id
1685 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 for(size_t i = 0; i < vUserData.size(); ++i)
1693 {
1694 for(size_t f = 0; f < TUserData::numFct; ++f)
1695 {
1696 // get function index
1697 const size_t fct = vUserData[i]->fct[f];
1698
1699 // get lfeID for function
1700 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 numSideIPs = err_est_data->get(fct)->num_side_ips(roid);
1710 sideLocIPs = err_est_data->get(fct)->template side_local_ips<side_type::dim>(roid);
1711 sideGlobIPs = err_est_data->get(fct)->side_global_ips(elem, &vCoCo[0]);
1712 }
1713 UG_CATCH_THROW("Global integration points for error estimator cannot be determined.");
1714
1715 // loop IPs
1716 for (size_t ip = 0; ip < numSideIPs; ++ip)
1717 {
1718 // get Dirichlet value (and do nothing, if conditional D value is false)
1719 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 (*err_est_data->get(fct))(elem,ip) = 0;
1725 continue;
1726 }
1727
1728 // evaluate shapes at ip
1729 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 rTrialSpace.shapes(vShape, sideLocIPs[ip]);
1734
1735 // get multiindices of element
1736 std::vector<DoFIndex> ind;
1737 dd->dof_indices(elem, fct, ind);
1738
1739 // compute solution at integration point
1740 number uAtIP = 0.0;
1741 for (size_t sh = 0; sh < vShape.size(); ++sh)
1742 uAtIP += DoFRef(u, ind[sh]) * vShape[sh];
1743
1744 // set error integrand value at IP
1745 (*err_est_data->get(fct))(elem,ip) = val[f] - uAtIP;
1746 }
1747 }
1748 }
1749 }
1750 }
1751 UG_CATCH_THROW("DirichletBoundary::adjust_error:"
1752 " While calling 'adjust_error' for TUserData, aborting.");
1753 }
1754}
1755
1756
1757} // end namespace ug
1758
1759#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__LAGRANGE_DIRICHLET_BOUNDARY_IMPL__ */
location name
Definition checkpoint_util.lua:128
Definition smart_pointer.h:296
Definition smart_pointer.h:108
base class for approximation spaces without type of algebra or dof distribution
Definition approximation_space.h:279
virtual void adjust_error(const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol=SPNULL, const std::vector< number > *vScaleMass=NULL, const std::vector< number > *vScaleStiff=NULL)
adjusts the assembled error estimator values in the attachments according to the constraint
Definition lagrange_dirichlet_boundary_impl.h:1599
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition lagrange_dirichlet_boundary.h:74
void add(SmartPtr< UserData< number, dim, bool > > func, const char *function, const char *subsets)
adds a lua callback (cond and non-cond)
Definition lagrange_dirichlet_boundary_impl.h:73
domain_type::position_type position_type
Type of position coordinates (e.g. position_type)
Definition lagrange_dirichlet_boundary.h:68
void assemble_dirichlet_rows(matrix_type &mat, ConstSmartPtr< DoFDistribution > dd, number time=0.0)
Sets dirichlet rows for all registered dirichlet values.
Definition lagrange_dirichlet_boundary_impl.h:430
algebra_type::vector_type vector_type
Type of algebra vector.
Definition lagrange_dirichlet_boundary.h:77
void adjust_linear(matrix_type &A, vector_type &b, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets unity rows in A and dirichlet values in right-hand side b
Definition lagrange_dirichlet_boundary_impl.h:1275
void adjust_defect(vector_type &d, const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol=NULL, const std::vector< number > *vScaleMass=NULL, const std::vector< number > *vScaleStiff=NULL)
sets a zero value in the defect for all dirichlet indices
Definition lagrange_dirichlet_boundary_impl.h:920
void set_approximation_space(SmartPtr< ApproximationSpace< TDomain > > approxSpace)
sets the approximation space to work on
Definition lagrange_dirichlet_boundary_impl.h:53
virtual void adjust_restriction(matrix_type &R, ConstSmartPtr< DoFDistribution > ddCoarse, ConstSmartPtr< DoFDistribution > ddFine, int type, number time=0.0)
sets constraints in restriction
Definition lagrange_dirichlet_boundary_impl.h:618
void check_functions_and_subsets(FunctionGroup &functionGroup, SubsetGroup &subsetGroup, size_t numFct) const
Definition lagrange_dirichlet_boundary_impl.h:309
MultipleSideAndElemErrEstData< TDomain > err_est_type
error estimator type
Definition lagrange_dirichlet_boundary.h:83
void adjust_rhs(vector_type &b, const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets the dirichlet value in the right-hand side
Definition lagrange_dirichlet_boundary_impl.h:1458
void adjust_solution(vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets the dirichlet value in the solution for all dirichlet indices
Definition lagrange_dirichlet_boundary_impl.h:1042
void adjust_jacobian(matrix_type &J, const vector_type &u, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0, ConstSmartPtr< VectorTimeSeries< vector_type > > vSol=NULL, const number s_a0=1.0)
sets a unity row for all dirichlet indices
Definition lagrange_dirichlet_boundary_impl.h:754
virtual void adjust_prolongation(matrix_type &P, ConstSmartPtr< DoFDistribution > ddFine, ConstSmartPtr< DoFDistribution > ddCoarse, int type, number time=0.0)
sets constraints in prolongation
Definition lagrange_dirichlet_boundary_impl.h:486
void clear()
removes all scheduled dirichlet data.
Definition lagrange_dirichlet_boundary_impl.h:63
void extract_data()
Definition lagrange_dirichlet_boundary_impl.h:410
virtual void adjust_correction(vector_type &c, ConstSmartPtr< DoFDistribution > dd, int type, number time=0.0)
sets zero to correction
Definition lagrange_dirichlet_boundary_impl.h:1159
Definition function_group.h:52
size_t size() const
number of functions in this group
Definition function_group.h:116
The base class for all geometric objects, such as vertices, edges, faces, volumes,...
Definition grid_base_objects.h:157
Identifier for Local Finite Elements.
Definition local_finite_element_id.h:98
SpaceType type() const
returns the type of the local finite element
Definition local_finite_element_id.h:132
int order() const
returns the order of the local finite element
Definition local_finite_element_id.h:126
int dim() const
returns the space dimension of the local finite element
Definition local_finite_element_id.h:129
virtual base class for local shape function sets
Definition local_shape_function_set.h:70
virtual void shapes(shape_type *vShape, const MathVector< dim > &x) const =0
returns all shape functions evaluated at a point
static SmartPtr< LuaUserData< TData, dim, TRet > > create(const std::string &name)
Definition lua_user_data.h:223
provides data specified in the lua script
Definition lua_user_data.h:96
static bool check_callback_returns(const char *callName, const bool bThrow=false)
returns true if callback has correct return values
Definition lua_user_data_impl.h:249
static std::string signature()
returns string of required callback signature
Definition lua_user_data_impl.h:67
a mathematical Vector with N entries.
Definition math_vector.h:97
Group of subsets.
Definition subset_group.h:51
size_t size() const
number of subsets in this group
Definition subset_group.h:122
Type based UserData.
Definition user_data.h:143
time series of solutions and corresponding time point
Definition solution_time_series.h:59
Base-class for all vertex-types.
Definition grid_base_objects.h:231
void SetRow(TSparseMatrix &A, size_t i, size_t alpha, number val=0.0)
Definition sparsematrix_util.h:710
void CollectCornerCoordinates(std::vector< typename TAAPos::ValueType > &vCornerCoordsOut, const TElem &elem, const TAAPos &aaPos, bool clearContainer=true)
returns the corner coordinates of a geometric object
Definition domain_util_impl.h:75
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
double number
Definition types.h:124
function util LuaCallbackHelper create(func)
the ug namespace
IndexLayout::Interface::iterator find(IndexLayout::Interface &interface, size_t i)
Definition parallel_index_layout.h:77
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition grid_base_objects.h:74
@ VOLUME
Definition grid_base_objects.h:63
@ VERTEX
Definition grid_base_objects.h:60
@ EDGE
Definition grid_base_objects.h:61
@ FACE
Definition grid_base_objects.h:62
bool CheckLuaCallbackName(const char *name)
returns true if callback exists
Definition lua_user_data.cpp:50
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition multi_index.h:276
function func(x, y, z, t, si)
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
grouping for subset and conditional data
Definition lagrange_dirichlet_boundary.h:363
grouping for subset and constant data
Definition lagrange_dirichlet_boundary.h:387
grouping for subset and non-conditional data
Definition lagrange_dirichlet_boundary.h:338
grouping for subset and the data already stored in the solution
Definition lagrange_dirichlet_boundary.h:435
grouping for subset and non-conditional vector data
Definition lagrange_dirichlet_boundary.h:411
SurfaceView::traits< TElem >::const_iterator const_iterator
Definition dof_distribution.h:82