ug4
active_set_impl.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3  * Author: Raphael Prohl
4  *
5  * This file is part of UG4.
6  *
7  * UG4 is free software: you can redistribute it and/or modify it under the
8  * terms of the GNU Lesser General Public License version 3 (as published by the
9  * Free Software Foundation) with the following additional attribution
10  * requirements (according to LGPL/GPL v3 §7):
11  *
12  * (1) The following notice must be displayed in the Appropriate Legal Notices
13  * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14  *
15  * (2) The following notice must be displayed at a prominent place in the
16  * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17  *
18  * (3) The following bibliography is recommended for citation and must be
19  * preserved in all covered files:
20  * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21  * parallel geometric multigrid solver on hierarchically distributed grids.
22  * Computing and visualization in science 16, 4 (2013), 151-164"
23  * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24  * flexible software system for simulating pde based models on high performance
25  * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26  *
27  * This program is distributed in the hope that it will be useful,
28  * but WITHOUT ANY WARRANTY; without even the implied warranty of
29  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  * GNU Lesser General Public License for more details.
31  */
32 
33 #ifndef ACTIVE_SET_IMPL_H_
34 #define ACTIVE_SET_IMPL_H_
35 
36 #include "active_set.h"
39 
40 using namespace std;
41 
42 namespace ug {
43 
44 template <int dim> struct face_type_traits
45 {
46  typedef void face_type0;
47  typedef void face_type1;
48  typedef void DimFEGeo;
49 };
50 
51 template <> struct face_type_traits<1>
52 {
56 };
57 
58 template <> struct face_type_traits<2>
59 {
63 };
64 
65 template <> struct face_type_traits<3>
66 {
70 };
71 
72 template <typename TDomain, typename TAlgebra>
74 {
75  m_vActiveSetGlob.resize(0); m_vActiveSetGlobOld.resize(0);
76  m_spDD = u.dof_distribution();
77  m_spDom = u.domain();
78 }
79 
80 /*template <typename TDomain, typename TAlgebra>
81 bool ActiveSet<TDomain, TAlgebra>::check_dist_to_obs(vector_type& u)
82 {
83  // STILL IN PROGRESS: u sollte hier reference-position + Startl�sung sein!
84  value_type dist;
85 
86  bool geometry_cut_by_cons = false;
87 
88  for(size_t i = 0; i < u.size(); i++)
89  {
90  UG_LOG("u(" << i << "):" << u[i] << "\n");
91  UG_LOG("m_spObs(" << i << "):" << (*m_spObs)[i] << "\n");
92  dist = (*m_spObs)[i] - u[i];
93  UG_LOG("dist:" << dist << "\n");
94  //TODO: anstatt u muss hier die geometrische Info einflie�en!
95  for (size_t fct = 0; fct < m_nrFcts; fct++)
96  {
97  if (BlockRef(dist,fct) < 0.0) // i.e.: u < m_spObs
98  {
99  geometry_cut_by_cons = true;
100  break;
101  }
102  }
103 
104  if (geometry_cut_by_cons)
105  break;
106 
107  }
108 
109  return geometry_cut_by_cons;
110 }*/
111 
112 // builds up a vector of global (dof,fct)-pairs, which are active (wrt. an obstacle-constraint)
113 template <typename TDomain, typename TAlgebra>
114 template <typename TElem, typename TIterator>
116  TIterator iterEnd,
117  function_type& u,
118  function_type& rhs,
119  function_type& lagrangeMult)
120 {
121  // check if at least an element exists, else return
122  if(iterBegin == iterEnd) return;
123 
124  int elemCounter = 0;
125 
126  static const int dim = function_type::dim;
127  typedef typename function_type::domain_type domain_type;
128  typedef typename face_type_traits<dim>::face_type0 face_type0;
129  typedef typename face_type_traits<dim>::face_type1 face_type1;
130 
131  // get position accessor
132  typename domain_type::position_accessor_type& aaPos
133  = u.domain()->position_accessor();
134 
135  // storage for corner coordinates
136  vector<MathVector<dim> > vCorner;
137  vector<MathVector<dim> > vSideCoPos;
138 
139  MathVector<dim> normal;
140 
141  // local indices and local algebra
142  LocalIndices ind, indObs;
143  LocalVector locU, locLM, locObs;
144  number complementaryVal;
145 
146  // TODO: it could be more efficient to store the active elements
147  // and its active local indices as well
148 
149  // Loop over all elements on active subsets
150  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
151  {
152  TElem* sideElem = *iter;
153 
154  // reference object type
155  ReferenceObjectID roid = sideElem->reference_object_id();
156 
157  const DimReferenceElement<dim-1>& rRefElem
159 
160  // get corners of element
161  CollectCornerCoordinates(vCorner, *sideElem, aaPos);
162 
163  // here the ordering of the corners in the reference element is exploited
164  // in order to compute the outer normal later on
165  int nCorner = (int)vCorner.size();
166  for (int i = 0; i < nCorner; ++i)
167  vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
168 
169  if (nCorner == dim)
170  {
171  ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
172  UG_LOG("face_type0 \n");
173  }
174  else
175  {
176  ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
177  UG_LOG("face_type1 \n");
178  }
179 
180  /*for (int i = 0; i < (int)vCorner.size(); ++i)
181  UG_LOG("sideCoPos: " << sideCoPos[i] << "\n");
182  UG_LOG("normal: " << normal << "\n");*/
183 
184  elemCounter++;
185 
186  // get global indices
187  u.indices(*iter, ind); (*m_spObs).indices(*iter, indObs);
188 
189  // adapt local algebra
190  locU.resize(ind); locLM.resize(ind); locObs.resize(indObs);
191 
192  // read local values of u and lagrangeMult
193  GetLocalVector(locU, u); GetLocalVector(locLM, lagrangeMult);
194  GetLocalVector(locObs, *m_spObs);
195 
196  // loop over DoFs in element and store all activeMultiIndex-pairs in vector
197  size_t nrFctElem = ind.num_fct();
198  number normOfNormal = VecLength(normal);
199  MathVector<dim> locUDof;
200 
201  for(size_t fct = 0; fct < nrFctElem; ++fct)
202  {
203  size_t nrDoFsPerFctElem = ind.num_dof(fct);
204  for(size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
205  {
206  for(int i = 0; i < dim; ++i)
207  locUDof[i] = locU(i, dof);
208 
209  // locLM corresponds to the lagrange multiplier lambda,
210  // c.f. Hintermueller/Ito/Kunisch(2003)
211  number locUNormal = VecDot(locUDof, normal) / normOfNormal;
212  complementaryVal = locLM(fct ,dof) + locUNormal - locObs(fct, dof);
213 
214  if (complementaryVal <= 1e-10){
215  locLM(fct,dof) = 0.0;
216  }
217  else
218  {
219  // create list of active global MultiIndex-pairs. Only those pairs should be attached
220  // which are not already a member of the 'activeSetGlob'-vector
221  bool bAlreadyActive = false;
222  for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
223  itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
224  {
225  if ((*itActiveInd)[0] == ind.index(fct, dof)
226  && (*itActiveInd)[1] == ind.comp(fct, dof))
227  bAlreadyActive = true;
228  }
229 
230  if (!bAlreadyActive)
231  {
232  DoFIndex newActiveIndex(ind.index(fct, dof), ind.comp(fct, dof));
233  m_vActiveSetGlob.push_back(newActiveIndex);
234 
235  BlockRef(rhs[ind.index(fct, dof)], ind.comp(fct, dof)) = locObs.value(fct,dof);
236  }
237  }
238  } // end(dof)
239  } // end(fct)
240 
241  // send local to global lagrangeMult
242  //AddLocalVector(lagrangeMult, locLM);
243 
244  UG_LOG("#activeDoFFctPairs global: " << m_vActiveSetGlob.size() << "\n");
245  } // end(elem)
246  UG_LOG("#elems: " << elemCounter << "\n");
247 
248 }
249 
250 // builds up a vector of global (dof,fct)-pairs, which are active (wrt. an obstacle-constraint)
251 // and sets dirichlet values in rhs for active indices
252 template <typename TDomain, typename TAlgebra>
254  function_type& rhs, function_type& lagrangeMult, function_type& gap)
255 {
256  // note: the active-index search is restricted
257  // to constraint of the form u * n <= consGF
258  if(!m_bObs)
259  UG_THROW("No constraint set in ActiveSet \n");
260 
261  if (u.num_indices() != rhs.num_indices() || u.num_indices() != lagrangeMult.num_indices() )
262  UG_THROW("GridFunctions u, rhs and lagrangeMult need to be defined "
263  "for the same domain and of the same size in 'ActiveSet:active_index' \n");
264 
265  // remember old ActiveSet for convergence check
266  // TODO: avoid this vector copy; m_vActiveSetGlobOld really necessary?
267  m_vActiveSetGlobOld = m_vActiveSetGlob;
268  m_vActiveSetGlob.resize(0);
269 
270  // 1.) get all subsets on which the 'gap'-gridfunction is defined!
271  // -> store them in vSubsetsOflagrangeMults
272  m_vActiveSubsets.resize(0);
273  //TODO: it is only necessary to loop over all BoundarySubsets!
274  for (int si = 0; si < m_spDD->num_subsets(); si++){
275  for (size_t fct = 0; fct < gap.num_fct(); fct++)
276  if (gap.is_def_in_subset(fct, si))
277  {
278  m_vActiveSubsets.push_back(si);
279  // 'break' is necessary to ensure that 'si' is
280  // only added once when several fcts of
281  // 'gap' are defined in subset 'si'!
282  break;
283  }
284  }
285 
286  if (m_vActiveSubsets.size() == 0)
287  UG_LOG("No subsets chosen as possible active subsets. \n");
288 
289  UG_LOG("#sizeOfvActiveSubsets: " << m_vActiveSubsets.size() << "\n");
290 
291  // 2.) loop over all elements of the possible active subsets
292  for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
293  activeSI != m_vActiveSubsets.end(); ++activeSI)
294  {
295  //UG_LOG("activeSI: " << *activeSI << "\n");
296  const int subsetDim = DimensionOfSubset(*m_spDD->subset_handler(), *activeSI);
297  //UG_LOG("subsetDim: " << subsetDim << "\n");
298 
299  // 3.) get localU out of u for each element and
300  // 4.) store the active global (dof,fct)-pair in m_vActiveSetGlob
301  switch(subsetDim)
302  {
303  case 0:
304  break;
305  case 1:
306  active_index_elem<RegularEdge>
307  (m_spDD->template begin<RegularEdge>(*activeSI),
308  m_spDD->template end<RegularEdge>(*activeSI), u, rhs, lagrangeMult);
309  break;
310  case 2:
311  active_index_elem<Triangle>
312  (m_spDD->template begin<Triangle>(*activeSI),
313  m_spDD->template end<Triangle>(*activeSI), u, rhs, lagrangeMult);
314  active_index_elem<Quadrilateral>
315  (m_spDD->template begin<Quadrilateral>(*activeSI),
316  m_spDD->template end<Quadrilateral>(*activeSI), u, rhs, lagrangeMult);
317  break;
318  default:
319  UG_THROW("ActiveSet::active_index:"
320  "SubsetDimension "<< subsetDim <<" (subset="<< *activeSI <<") not supported.");
321  }
322  }
323 
324  if (m_vActiveSetGlob.size() > 0) return true;
325  else return false;
326 }
327 
328 // sets a Dirichlet row for all active Indices
329 template <typename TDomain, typename TAlgebra>
332 {
333  for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
334  itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd){
335  SetDirichletRow(mat, *itActiveInd);
336  }
337 }
338 
339 // computes the lagrange multiplier for a given discretization
340 template <typename TDomain, typename TAlgebra>
342  const function_type& u)
343 {
344  // check that lagrange multiplier disc is set
345  if (m_spLagMultDisc.invalid())
346  UG_THROW("No discretization set to compute the lagrange multiplier (in "
347  "ActiveSet:lagrange_multiplier) \n");
348 
349  if(m_vActiveSetGlob.size() != 0.0)
350  {
351  UG_LOG("activeDoFs in ActiveSet:lagrange_multiplier " << m_vActiveSetGlob.size() << "\n");
352  //TODO: pass localActiveElemAndIndex here!!!
353  m_spLagMultDisc->lagrange_multiplier(lagrangeMult, u, m_vActiveSetGlob, m_vActiveSubsets);
354  }
355 }
356 
357 template <typename TDomain, typename TAlgebra>
358 template <typename TElem, typename TIterator>
360  TIterator iterEnd, matrix_type& lagrangeMatInv)
361 {
362  typedef typename face_type_traits<dim>::face_type0 face_type0;
363  typedef typename face_type_traits<dim>::face_type1 face_type1;
364  typedef typename face_type_traits<dim>::DimFEGeo sideGeo;
365  typename TDomain::position_accessor_type& aaPos = m_spDom->position_accessor();
366 
367  // some storage
368  MathVector<dim> normal;
369  vector<MathVector<dim> > vCorner;
370  vector<MathVector<dim> > vSideCoPos;
371 
372  // local indices and local algebra
373  LocalIndices ind; LocalMatrix loclagrangeMatInv;
374 
375  // Loop over all elements on active subsets
376  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
377  {
378  TElem* elem = *iter;
379 
380  // get global indices
381  m_spDD->indices(elem, ind);
382 
383  loclagrangeMatInv.resize(ind, ind);
384  loclagrangeMatInv = 0.0;
385 
386  // reference object type and geometry
387  ReferenceObjectID sideRoid = elem->reference_object_id();
388  sideGeo geo(sideRoid, 3, LFEID(LFEID::LAGRANGE, dim, 1));
389 
390  // prepare geometry for type and order
391  try{
392  geo.update_local(sideRoid, LFEID(LFEID::LAGRANGE, dim, 1), 1);
393  }
394  UG_CATCH_THROW("ActiveSet::lagrange_mat_inv_elem:"
395  " Cannot update local values of finite element geometry.");
396 
397  // get corners of element
398  CollectCornerCoordinates(vCorner, *elem, aaPos);
399 
400  const DimReferenceElement<dim-1>& rRefElem
401  = ReferenceElementProvider::get<dim-1>(sideRoid);
402 
403  int nCorner = (int)vCorner.size();
404  for (int i = 0; i < nCorner; ++i)
405  vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
406 
407  if (nCorner == dim)
408  ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
409  else
410  ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
411 
412  UG_LOG("normal in lagrange_mat_inv_elem: " << normal << "\n");
413  //number normOfNormal = VecLength(normal);
414 
415  try{
416  geo.update(elem, &vSideCoPos[0], LFEID(LFEID::LAGRANGE, dim, 1), 1);
417  }
418  UG_CATCH_THROW("ActiveSet::lagrange_mat_inv_elem:"
419  " Cannot update finite element geometry.");
420 
421  for (size_t ip = 0; ip < geo.num_ip(); ++ip)
422  for(size_t sh = 0; sh < geo.num_sh(); ++sh)
423  for (size_t i = 0; i < (size_t) dim; ++i)
424  loclagrangeMatInv(i, sh, i, sh) += 1.0/(geo.weight(ip) * geo.shape(ip, sh));
425  // * normal[i] / normOfNormal;
426 
427  AddLocalMatrixToGlobal(lagrangeMatInv, loclagrangeMatInv);
428  }
429 }
430 
431 template <typename TDomain, typename TAlgebra>
433 {
434  for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
435  activeSI != m_vActiveSubsets.end(); ++activeSI)
436  {
437  UG_LOG("activeSI: " << *activeSI << "\n");
438  const int subsetDim = DimensionOfSubset(*m_spDD->subset_handler(), *activeSI);
439  UG_LOG("subsetDim: " << subsetDim << "\n");
440 
441  switch(subsetDim)
442  {
443  case 0:
444  break;
445  case 1:
446  lagrange_mat_inv_elem<RegularEdge>
447  (m_spDD->template begin<RegularEdge>(*activeSI), m_spDD->template end<RegularEdge>(*activeSI), lagrangeMatInv);
448  break;
449  case 2:
450  lagrange_mat_inv_elem<Triangle>
451  (m_spDD->template begin<Triangle>(*activeSI), m_spDD->template end<Triangle>(*activeSI), lagrangeMatInv);
452  lagrange_mat_inv_elem<Quadrilateral>
453  (m_spDD->template begin<Quadrilateral>(*activeSI), m_spDD->template end<Quadrilateral>(*activeSI), lagrangeMatInv);
454  break;
455  default:
456  UG_THROW("ActiveSet::lagrange_mat_inv:"
457  "SubsetDimension "<<subsetDim<<" (subset="<<*activeSI<<") not supported.");
458  }
459  }
460 }
461 
462 template <typename TDomain, typename TAlgebra>
464  const matrix_type& mat,
465  const vector_type& u,
466  vector_type& rhs)
467 {
468  // only if some indices are active the lagrange multiplier is computed
469  if(m_vActiveSetGlob.size() != 0.0)
470  {
471  if (u.size() != lagMult.size())
472  UG_THROW("Temporarily u and lagMult need to be "
473  "of same size in ActiveSet:residual_lagrange_mult \n");
474 
475  /*matrix_type lagrangeMatInv;
476  //SmartPtr<AssembledLinearOperator<algebra_type> > splagrangeMatInv;
477  ass_lagrangeMatInv(lagrangeMatInv);*/
478 
479  SmartPtr<vector_type> spMat_u = u.clone_without_values();
480  (*spMat_u).resize(u.size());
481 
482  #ifdef UG_PARALLEL
483  MatMultDirect(*spMat_u, 1.0, mat, u);
484  spMat_u->set_storage_type(u.get_storage_mask());
485  #else
486  MatMult(*spMat_u, 1.0, mat, u);
487  #endif
488 
489  //SmartPtr<vector_type> spRes = u.clone_without_values();
490  //(*spRes).resize(u.size());
491 
492  // loop MultiIndex-pairs in activeSet-vector
493  for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
494  itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
495  {
496  // compute lagrange multiplier for active multiIndices
497  // lagMult = rhs - Mat * u;
498  DoFRef(lagMult, *itActiveInd) = DoFRef(rhs, *itActiveInd) - DoFRef((*spMat_u), *itActiveInd);
499  }
500 
501  /*#ifdef UG_PARALLEL
502  MatMultDirect(lagMult, 1.0, lagrangeMatInv, *spRes);
503  #else
504  MatMult(lagMult, 1.0, lagrangeMatInv, *spRes);
505  #endif*/
506 
507  UG_LOG("new lagMult-values computed \n");
508  //UG_LOG("rhs updated \n");
509  }
510  else{
511  UG_LOG("no active index in residual_lagrange_mult \n");
512  }
513 }
514 
515 
516 template <typename TDomain, typename TAlgebra>
517 template <typename TElem, typename TIterator>
519  TIterator iterEnd, function_type& u, const function_type& lambda)
520 {
521  static const int dim = function_type::dim;
522  typedef typename function_type::domain_type domain_type;
523  typedef typename face_type_traits<dim>::face_type0 face_type0;
524  typedef typename face_type_traits<dim>::face_type1 face_type1;
525 
526  // get position accessor
527  typename domain_type::position_accessor_type& aaPos
528  = u.domain()->position_accessor();
529 
530  // storage for corner coordinates
531  vector<MathVector<dim> > vCorner;
532  vector<MathVector<dim> > vSideCoPos;
533  MathVector<dim> normal;
534 
535  // local indices and local algebra
536  LocalIndices ind, indObs;
537  LocalVector locU, locObs, locLambda;
538 
539  // Loop over all elements on active subsets
540  for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
541  {
542  TElem* elem = *iter;
543 
544  // reference object type
545  ReferenceObjectID roid = elem->reference_object_id();
546 
547  const DimReferenceElement<dim-1>& rRefElem
549 
550  // get corners of element
551  CollectCornerCoordinates(vCorner, *elem, aaPos);
552 
553  // here the ordering of the corners in the reference element is exploited
554  // in order to compute the outer normal later on
555  int nCorner = (int)vCorner.size();
556  for (int i = 0; i < nCorner; ++i)
557  vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
558 
559  if (nCorner == dim)
560  ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
561  else
562  ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
563 
564  // get global indices
565  u.indices(*iter, ind); (*m_spObs).indices(*iter, indObs);
566 
567  // adapt local algebra
568  locU.resize(ind); locObs.resize(indObs); locLambda.resize(ind);
569 
570  // read local values of u and lagrangeMult
571  GetLocalVector(locU, u); GetLocalVector(locObs, *m_spObs);
572  GetLocalVector(locLambda, lambda);
573 
574  size_t nrFctElem = ind.num_fct();
575  number gapValue, locUNormal, kktcond;
576  number normOfNormal = VecLength(normal);
577  MathVector<dim> locUDof;
578 
579  for(size_t fct = 0; fct < nrFctElem; ++fct)
580  {
581  size_t nrDoFsPerFctElem = ind.num_dof(fct);
582  for(size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
583  {
584  for(int i = 0; i < dim; ++i)
585  locUDof[i] = locU(i, dof);
586  locUNormal = VecDot(locUDof, normal) / normOfNormal;
587 
588  gapValue = locUNormal - locObs(fct, dof);
589  if (gapValue > 1e-06){
590  // i.e.: m_spObs < u
591  // constraint is violated
592  return false;
593  }
594 
595  kktcond = gapValue * locLambda(fct, dof);
596  if (kktcond > 1e-06 || kktcond < -1e-06)
597  return false;
598 
599  } // end(dof)
600  } // end(fct)
601  } // end(elem)
602 
603  return true;
604 }
605 
606 template <typename TDomain, typename TAlgebra>
608  const size_t step)
609 {
610  // ensure that at least one activeSet-iteration is performed
611  if (step <= 1)
612  return false;
613 
614  // NOW TWO CHECKS WILL BE PERFORMED TO ENSURE CONVERGENCE:
615  // 1. Did some multiIndices change from 'active' to 'inactive' or vice versa
616  // in the last iteration-step?
617  // 2. Is the constraint violated for any DoFIndex?
618 
619  UG_LOG(m_vActiveSetGlob.size() << " indices are active at the begin "
620  "of step " << step << " ! \n");
621 
622  // check if activeSet has changed
623  if (m_vActiveSetGlob == m_vActiveSetGlobOld)
624  {
625  // check if constraint is fulfilled
626  bool bConstraintViolated = false;
627  for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
628  activeSI != m_vActiveSubsets.end(); ++activeSI)
629  {
630  const int subsetDim = DimensionOfSubset(*m_spDD->subset_handler(), *activeSI);
631  switch(subsetDim)
632  {
633  case 0:
634  break;
635  case 1:
636  if (!check_conv_elem<RegularEdge>
637  (m_spDD->template begin<RegularEdge>(*activeSI), m_spDD->template end<RegularEdge>(*activeSI),
638  u, lambda))
639  {bConstraintViolated = true;}
640 
641  break;
642  case 2:
643  if (!check_conv_elem<Triangle>
644  (m_spDD->template begin<Triangle>(*activeSI),
645  m_spDD->template end<Triangle>(*activeSI), u, lambda))
646  {bConstraintViolated = true;}
647 
648  if (!check_conv_elem<Quadrilateral>
649  (m_spDD->template begin<Quadrilateral>(*activeSI),
650  m_spDD->template end<Quadrilateral>(*activeSI), u, lambda))
651  {bConstraintViolated = true;}
652 
653  break;
654  default:
655  UG_THROW("ActiveSet::check_conv:"
656  "SubsetDimension "<< subsetDim <<" (subset="<< *activeSI <<") not supported.");
657  }
658 
659  if (bConstraintViolated)
660  return false;
661  }
662 
663  // activeSet remains unchanged & constraint is fulfilled for all indices
664  return true;
665  }
666  else{
667  return false;
668  }
669 }
670 
671 
672 template <typename TDomain, typename TAlgebra>
674  const vector_type& u,
675  const vector_type& lambda,
676  const vector_type& rhs)
677 {
678  if (u.size() != lambda.size())
679  UG_THROW("Temporarily u and lambda need to be "
680  "of same size in ActiveSet:check_inequ \n");
681 
682  SmartPtr<vector_type> spMat_u = u.clone_without_values();
683  SmartPtr<vector_type> spRes = u.clone_without_values();
684  (*spMat_u).resize(u.size()); (*spRes).resize(u.size());
685 
686  #ifdef UG_PARALLEL
687  MatMultDirect(*spMat_u, 1.0, mat, u);
688  spMat_u->set_storage_type(u.get_storage_mask());
689  spRes->set_storage_type(u.get_storage_mask());
690  #else
691  MatMult(*spMat_u, 1.0, mat, u);
692  #endif
693 
694  for (size_t i = 0; i < u.size(); ++i)
695  {
696  //if (lambda[i] < -1e-06)
697  // return false;
698 
699  (*spRes)[i] = (*spMat_u)[i] - rhs[i]; //+ lambda[i] - rhs[i]; //
700  UG_LOG("lambda["<< i << "]: " << lambda[i] <<", res[" << i << "]: " << (*spRes)[i] << "\n");
701  //if ((*spRes)[i] > 1e-06)
702  // return false;
703  }
704 
705  return true;
706 }
707 
708 }; // namespace ug
709 
710 #endif /* ACTIVE_SET_IMPL_H_ */
bool MatMult(vector_t &dest, const matrix_t &A1, const vector_t &w1)
calculates dest = A1 * w1;
Active Set method.
Definition: active_set.h:100
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition: active_set.h:106
algebra_type::vector_type vector_type
Type of algebra vector.
Definition: active_set.h:109
Definition: fe_geom.h:183
dimension dependent base class for reference elements
Definition: reference_element.h:183
bool is_def_in_subset(size_t fct, int si) const
returns if a function is defined on a subset
Definition: dof_distribution_info.h:217
size_t num_fct() const
number of discrete functions on subset si
Definition: dof_distribution_info.h:176
represents numerical solutions on a grid using an algebraic vector
Definition: grid_function.h:121
TDomain domain_type
Domain.
Definition: grid_function.h:130
size_t num_indices() const
return the number of indices distributed (proc local)
Definition: grid_function.h:303
SmartPtr< TDomain > domain()
returns domain
Definition: grid_function.h:342
void indices(TElem *elem, LocalIndices &ind, bool bHang=false) const
get all indices of the element
Definition: grid_function.h:317
SmartPtr< DoFDistribution > dof_distribution()
Definition: grid_function.h:236
Identifier for Local Finite Elements.
Definition: local_finite_element_id.h:98
Definition: local_algebra.h:50
size_t num_dof(size_t fct) const
number of dofs for accessible function
Definition: local_algebra.h:117
index_type index(size_t fct, size_t dof) const
global algebra index for (fct, dof)
Definition: local_algebra.h:139
size_t num_fct() const
number of functions
Definition: local_algebra.h:114
comp_type comp(size_t fct, size_t dof) const
algebra comp for (fct, dof)
Definition: local_algebra.h:153
Definition: local_algebra.h:422
void resize(const LocalIndices &ind)
resize for same ind in row and column
Definition: local_algebra.h:445
Definition: local_algebra.h:198
void resize(const LocalIndices &ind)
resize for current local indices
Definition: local_algebra.h:214
Definition: multi_index.h:50
Definition: reference_element.h:345
Definition: reference_element.h:445
Definition: reference_element.h:394
Definition: reference_element.h:303
void SetDirichletRow(TSparseMatrix &A, size_t i, size_t alpha)
Definition: sparsematrix_util.h:796
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
#define UG_CATCH_THROW(msg)
Definition: error.h:64
#define UG_THROW(msg)
Definition: error.h:57
#define UG_LOG(msg)
Definition: log.h:367
double number
Definition: types.h:124
vector_t::value_type VecLength(const vector_t &v)
returns the length of v. Slower than VecLengthSq.
Definition: math_vector_functions_common_impl.hpp:341
vector_t::value_type VecDot(const vector_t &v1, const vector_t &v2)
returns the dot-product of two vector_ts
Definition: math_vector_functions_common_impl.hpp:385
SM_edge_weight_map< typename T::value_type, ug::BidirectionalMatrix< T > > get(edge_weight_t, ug::BidirectionalMatrix< T > const &g)
Definition: bidirectional_boost.h:157
Definition: smart_pointer.h:814
the ug namespace
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
double & BlockRef(T &vec, size_t i)
Definition: blocks.h:66
bool MatMultDirect(ParallelVector< vector_type > &dest, const number &beta1, const ParallelMatrix< matrix_type > &A1, const ParallelVector< vector_type > &w1)
Definition: parallel_matrix_impl.h:224
void GetLocalVector(LocalVector &lvec, const TVector &vec)
Definition: local_algebra.h:739
int DimensionOfSubset(const ISubsetHandler &sh, int si)
returns the current dimension of the subset
Definition: subset_dim_util.cpp:89
void AddLocalMatrixToGlobal(TMatrix &mat, const LocalMatrix &lmat)
Definition: local_algebra.h:767
ReferenceVertex face_type1
Definition: active_set_impl.h:54
ReferenceVertex face_type0
Definition: active_set_impl.h:53
DimFEGeometry< 1, 1 > DimFEGeo
Definition: active_set_impl.h:55
ReferenceEdge face_type0
Definition: active_set_impl.h:60
ReferenceEdge face_type1
Definition: active_set_impl.h:61
DimFEGeometry< 2, 1 > DimFEGeo
Definition: active_set_impl.h:62
ReferenceTriangle face_type0
Definition: active_set_impl.h:67
ReferenceQuadrilateral face_type1
Definition: active_set_impl.h:68
DimFEGeometry< 3, 2 > DimFEGeo
Definition: active_set_impl.h:69
Definition: active_set_impl.h:45
void face_type0
Definition: active_set_impl.h:46
void face_type1
Definition: active_set_impl.h:47
void DimFEGeo
Definition: active_set_impl.h:48