ug4
Loading...
Searching...
No Matches
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
40using namespace std;
41
42namespace ug {
43
44template <int dim> struct face_type_traits
45{
46 typedef void face_type0;
47 typedef void face_type1;
48 typedef void DimFEGeo;
49};
50
57
64
71
72template <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>
81bool 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)
113template <typename TDomain, typename TAlgebra>
114template <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
158 = ReferenceElementProvider::get<dim-1>(roid);
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
252template <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
329template <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
340template <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
357template <typename TDomain, typename TAlgebra>
358template <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
431template <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
462template <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
516template <typename TDomain, typename TAlgebra>
517template <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
548 = ReferenceElementProvider::get<dim-1>(roid);
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
606template <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
672template <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_ */
Definition smart_pointer.h:108
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition active_set.h:106
void lagrange_mat_inv(matrix_type &lagrangeMatInv)
Definition active_set_impl.h:432
bool active_index(function_type &u, function_type &rhs, function_type &lagrangeMult, function_type &gap)
determines the active indices, stores them in a vector and sets dirichlet values in rhs for active in...
Definition active_set_impl.h:253
bool check_conv(function_type &u, const function_type &lambda, const size_t step)
checks if all constraints are fulfilled & the activeSet remained unchanged
Definition active_set_impl.h:607
bool check_conv_elem(TIterator iterBegin, TIterator iterEnd, function_type &u, const function_type &lambda)
Definition active_set_impl.h:518
algebra_type::vector_type vector_type
Type of algebra vector.
Definition active_set.h:109
bool check_inequ(const matrix_type &mat, const vector_type &u, const vector_type &lagrangeMult, const vector_type &rhs)
checks if all inequalities are fulfilled
Definition active_set_impl.h:673
void residual_lagrange_mult(vector_type &lagMult, const matrix_type &mat, const vector_type &u, vector_type &rhs)
Definition active_set_impl.h:463
void lagrange_multiplier(function_type &lagrangeMult, const function_type &u)
computes the lagrange multiplier for a given disc
Definition active_set_impl.h:341
void lagrange_mat_inv_elem(TIterator iterBegin, TIterator iterEnd, matrix_type &lagrangeMatInv)
Definition active_set_impl.h:359
void active_index_elem(TIterator iterBegin, TIterator iterEnd, function_type &u, function_type &rhs, function_type &lagrangeMult)
checks the distance to the prescribed obstacle/constraint
Definition active_set_impl.h:115
void set_dirichlet_rows(matrix_type &mat)
Definition active_set_impl.h:331
void prepare(function_type &u)
Definition active_set_impl.h:73
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
void indices(TElem *elem, LocalIndices &ind, bool bHang=false) const
get all indices of the element
Definition grid_function.h:317
SmartPtr< TDomain > domain()
returns domain
Definition grid_function.h:342
SmartPtr< DoFDistribution > dof_distribution()
Definition grid_function.h:236
Identifier for Local Finite Elements.
Definition local_finite_element_id.h:98
@ LAGRANGE
Definition local_finite_element_id.h:104
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
a mathematical Vector with N entries.
Definition math_vector.h:97
Definition multi_index.h:50
Definition reference_element.h:345
number size() const
returns the size (e.g. area or volume) of the reference element
Definition reference_element.h:83
static const DimReferenceElement< dim > & get(ReferenceObjectID roid)
returns a dimension dependent Reference Element
Definition reference_element.h:280
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
#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
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
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
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
void AddLocalMatrixToGlobal(TMatrix &mat, const LocalMatrix &lmat)
Definition local_algebra.h:767
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition multi_index.h:276
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