ug4
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 
44 namespace ug{
45 
46 
48 // setup
50 
51 template <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 
61 template <typename TDomain, typename TAlgebra>
63 clear()
64 {
65  m_vBNDNumberData.clear();
66  m_vNumberData.clear();
67  m_vConstNumberData.clear();
68  m_vVectorData.clear();
69 }
70 
71 template <typename TDomain, typename TAlgebra>
73 add(SmartPtr<UserData<number, dim, bool> > func, const char* function, const char* subsets)
74 {
75  m_vBNDNumberData.push_back(CondNumberData(func, function, subsets));
76 }
77 
78 template <typename TDomain, typename 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  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 
96 template <typename TDomain, typename TAlgebra>
98 add(SmartPtr<UserData<number, dim> > func, const char* function, const char* subsets)
99 {
100  m_vNumberData.push_back(NumberData(func, function, subsets));
101 }
102 
103 template <typename TDomain, typename 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  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 
121 template <typename TDomain, typename TAlgebra>
123 add(number value, const char* function, const char* subsets)
124 {
125  m_vConstNumberData.push_back(ConstNumberData(value, function, subsets));
126 }
127 
128 template <typename TDomain, typename TAlgebra>
130 add(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 
146 template <typename TDomain, typename TAlgebra>
148 add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions, const char* subsets)
149 {
150  m_vVectorData.push_back(VectorData(func, functions, subsets));
151 }
152 
153 template <typename TDomain, typename 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  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
172 template <typename TDomain, typename TAlgebra>
174 add(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"
204  << (LuaUserData<number, dim>::signature()) << "\n" <<
205  "b) Conditional Number - Callback\n"
207  "c) "<<dim<<"d Vector - Callback\n"
208  << (LuaUserData<MathVector<dim>, dim>::signature()));
209 }
210 
211 template <typename TDomain, typename TAlgebra>
213 add(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 
229 template <typename TDomain, typename TAlgebra>
231 add(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"
256  << (LuaUserData<number, dim>::signature()) << "\n" <<
257  "b) Conditional Number - Callback\n"
259  "c) "<<dim<<"d Vector - Callback\n"
260  << (LuaUserData<MathVector<dim>, dim>::signature()));
261 }
262 
263 template <typename TDomain, typename TAlgebra>
265 add(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 
282 template <typename TDomain, typename TAlgebra>
284 add(const char* functions, const char* subsets)
285 {
286  m_vOldNumberData.push_back(OldNumberData(functions, subsets));
287 }
288 
289 template <typename TDomain, typename TAlgebra>
291 add(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 
307 template <typename TDomain, typename TAlgebra>
309 check_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 
353 template <typename TDomain, typename TAlgebra>
354 template <typename TUserData, typename TScheduledUserData>
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  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 
408 template <typename TDomain, typename TAlgebra>
410 extract_data()
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 
428 template <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 
441  DoFDistribution::traits<Vertex>::const_iterator iterBegin = dd->begin<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 
480 // adjust TRANSFER
482 
483 
484 template <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 
510 template <typename TDomain, typename TAlgebra>
511 template <typename TUserData>
513 adjust_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 
542 template <typename TDomain, typename TAlgebra>
543 template <typename TBaseElem, typename TUserData>
545 adjust_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 
616 template <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 
642 template <typename TDomain, typename TAlgebra>
643 template <typename TUserData>
645 adjust_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 
674 template <typename TDomain, typename TAlgebra>
675 template <typename TBaseElem, typename TUserData>
677 adjust_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 
752 template <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 
771 template <typename TDomain, typename TAlgebra>
772 template <typename TUserData>
774 adjust_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 
805 template <typename TDomain, typename TAlgebra>
806 template <typename TBaseElem, typename TUserData>
808 adjust_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 
918 template <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 
938 template <typename TDomain, typename TAlgebra>
939 template <typename TUserData>
941 adjust_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 
972 template <typename TDomain, typename TAlgebra>
973 template <typename TBaseElem, typename TUserData>
975 adjust_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 
1040 template <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 
1056 template <typename TDomain, typename TAlgebra>
1057 template <typename TUserData>
1059 adjust_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 
1089 template <typename TDomain, typename TAlgebra>
1090 template <typename TBaseElem, typename TUserData>
1092 adjust_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 
1157 template <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 
1177 template <typename TDomain, typename TAlgebra>
1178 template <typename TUserData>
1180 adjust_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 
1210 template <typename TDomain, typename TAlgebra>
1211 template <typename TBaseElem, typename TUserData>
1213 adjust_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 
1273 template <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 
1290 template <typename TDomain, typename TAlgebra>
1291 template <typename TUserData>
1293 adjust_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 
1324 template <typename TDomain, typename TAlgebra>
1325 template <typename TBaseElem, typename TUserData>
1327 adjust_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 
1456 template <typename TDomain, typename TAlgebra>
1458 adjust_rhs(vector_type& b, const vector_type& u,
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 
1473 template <typename TDomain, typename TAlgebra>
1474 template <typename TUserData>
1476 adjust_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 
1507 template <typename TDomain, typename TAlgebra>
1508 template <typename TBaseElem, typename TUserData>
1510 adjust_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 
1596 template <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 
1635 template <typename TDomain, typename TAlgebra>
1636 template <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
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
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
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
algebra_type::vector_type vector_type
Type of algebra vector.
Definition: constraint_interface.h:72
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition: constraint_interface.h:69
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:247
static std::string signature()
returns string of required callback signature
Definition: lua_user_data_impl.h:65
TErrEstData * get(size_t uid)
accessing the underlying error estimator data objects via function id
Definition: err_est_data.h:466
Definition: err_est_data.h:506
SideAndElemErrEstData< TDomain >::side_type side_type
type of the sides (face, edge) and the elems (volume, face)
Definition: err_est_data.h:512
Group of subsets.
Definition: subset_group.h:51
size_t size() const
number of subsets in this group
Definition: subset_group.h:122
Const iterator to traverse the surface of a multi-grid hierarchy.
Definition: surface_view.h:237
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
static const int dim
base_type::TBaseElem TBaseElem
#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
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition: multi_index.h:276
@ 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
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
T value_type
Definition: sparsematrix_interface.h:2
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