Line data Source code
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"
37 : #include "lib_disc/common/geometry_util.h"
38 : #include "lib_disc/spatial_disc/disc_util/fe_geom.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 : {
53 : typedef ReferenceVertex face_type0;
54 : typedef ReferenceVertex face_type1;
55 : typedef DimFEGeometry<1, 1> DimFEGeo;
56 : };
57 :
58 : template <> struct face_type_traits<2>
59 : {
60 : typedef ReferenceEdge face_type0;
61 : typedef ReferenceEdge face_type1;
62 : typedef DimFEGeometry<2, 1> DimFEGeo;
63 : };
64 :
65 : template <> struct face_type_traits<3>
66 : {
67 : typedef ReferenceTriangle face_type0;
68 : typedef ReferenceQuadrilateral face_type1;
69 : typedef DimFEGeometry<3, 2> DimFEGeo;
70 : };
71 :
72 : template <typename TDomain, typename TAlgebra>
73 0 : void ActiveSet<TDomain, TAlgebra>::prepare(function_type& u)
74 : {
75 0 : m_vActiveSetGlob.resize(0); m_vActiveSetGlobOld.resize(0);
76 0 : m_spDD = u.dof_distribution();
77 0 : m_spDom = u.domain();
78 0 : }
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>
115 0 : void ActiveSet<TDomain, TAlgebra>::active_index_elem(TIterator iterBegin,
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 0 : 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
151 : {
152 : TElem* sideElem = *iter;
153 :
154 : // reference object type
155 0 : 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 0 : int nCorner = (int)vCorner.size();
166 0 : for (int i = 0; i < nCorner; ++i)
167 0 : vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
168 :
169 0 : if (nCorner == dim)
170 : {
171 0 : ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
172 : UG_LOG("face_type0 \n");
173 : }
174 : else
175 : {
176 0 : 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 0 : elemCounter++;
185 :
186 : // get global indices
187 : u.indices(*iter, ind); (*m_spObs).indices(*iter, indObs);
188 :
189 : // adapt local algebra
190 0 : locU.resize(ind); locLM.resize(ind); locObs.resize(indObs);
191 :
192 : // read local values of u and lagrangeMult
193 0 : GetLocalVector(locU, u); GetLocalVector(locLM, lagrangeMult);
194 0 : 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 0 : for(size_t fct = 0; fct < nrFctElem; ++fct)
202 : {
203 : size_t nrDoFsPerFctElem = ind.num_dof(fct);
204 0 : for(size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
205 : {
206 0 : for(int i = 0; i < dim; ++i)
207 0 : locUDof[i] = locU(i, dof);
208 :
209 : // locLM corresponds to the lagrange multiplier lambda,
210 : // c.f. Hintermueller/Ito/Kunisch(2003)
211 0 : number locUNormal = VecDot(locUDof, normal) / normOfNormal;
212 0 : complementaryVal = locLM(fct ,dof) + locUNormal - locObs(fct, dof);
213 :
214 0 : if (complementaryVal <= 1e-10){
215 0 : 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 0 : for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
223 0 : itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
224 : {
225 0 : if ((*itActiveInd)[0] == ind.index(fct, dof)
226 0 : && (*itActiveInd)[1] == ind.comp(fct, dof))
227 : bAlreadyActive = true;
228 : }
229 :
230 0 : if (!bAlreadyActive)
231 : {
232 0 : DoFIndex newActiveIndex(ind.index(fct, dof), ind.comp(fct, dof));
233 0 : m_vActiveSetGlob.push_back(newActiveIndex);
234 :
235 0 : 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 0 : UG_LOG("#elems: " << elemCounter << "\n");
247 :
248 0 : }
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>
253 0 : bool ActiveSet<TDomain, TAlgebra>::active_index(function_type& u,
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 0 : if(!m_bObs)
259 0 : UG_THROW("No constraint set in ActiveSet \n");
260 :
261 0 : if (u.num_indices() != rhs.num_indices() || u.num_indices() != lagrangeMult.num_indices() )
262 0 : 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 0 : m_vActiveSetGlobOld = m_vActiveSetGlob;
268 0 : m_vActiveSetGlob.resize(0);
269 :
270 : // 1.) get all subsets on which the 'gap'-gridfunction is defined!
271 : // -> store them in vSubsetsOflagrangeMults
272 0 : m_vActiveSubsets.resize(0);
273 : //TODO: it is only necessary to loop over all BoundarySubsets!
274 0 : for (int si = 0; si < m_spDD->num_subsets(); si++){
275 0 : for (size_t fct = 0; fct < gap.num_fct(); fct++)
276 0 : if (gap.is_def_in_subset(fct, si))
277 : {
278 0 : 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 0 : 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 0 : for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
293 0 : activeSI != m_vActiveSubsets.end(); ++activeSI)
294 : {
295 : //UG_LOG("activeSI: " << *activeSI << "\n");
296 0 : 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 0 : switch(subsetDim)
302 : {
303 : case 0:
304 : break;
305 0 : case 1:
306 : active_index_elem<RegularEdge>
307 0 : (m_spDD->template begin<RegularEdge>(*activeSI),
308 : m_spDD->template end<RegularEdge>(*activeSI), u, rhs, lagrangeMult);
309 0 : break;
310 0 : case 2:
311 : active_index_elem<Triangle>
312 0 : (m_spDD->template begin<Triangle>(*activeSI),
313 : m_spDD->template end<Triangle>(*activeSI), u, rhs, lagrangeMult);
314 : active_index_elem<Quadrilateral>
315 0 : (m_spDD->template begin<Quadrilateral>(*activeSI),
316 : m_spDD->template end<Quadrilateral>(*activeSI), u, rhs, lagrangeMult);
317 0 : break;
318 0 : default:
319 0 : UG_THROW("ActiveSet::active_index:"
320 : "SubsetDimension "<< subsetDim <<" (subset="<< *activeSI <<") not supported.");
321 : }
322 : }
323 :
324 0 : if (m_vActiveSetGlob.size() > 0) return true;
325 0 : else return false;
326 : }
327 :
328 : // sets a Dirichlet row for all active Indices
329 : template <typename TDomain, typename TAlgebra>
330 0 : void ActiveSet<TDomain, TAlgebra>::
331 : set_dirichlet_rows(matrix_type& mat)
332 : {
333 0 : for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
334 0 : itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd){
335 : SetDirichletRow(mat, *itActiveInd);
336 : }
337 0 : }
338 :
339 : // computes the lagrange multiplier for a given discretization
340 : template <typename TDomain, typename TAlgebra>
341 0 : void ActiveSet<TDomain, TAlgebra>::lagrange_multiplier(function_type& lagrangeMult,
342 : const function_type& u)
343 : {
344 : // check that lagrange multiplier disc is set
345 0 : if (m_spLagMultDisc.invalid())
346 0 : UG_THROW("No discretization set to compute the lagrange multiplier (in "
347 : "ActiveSet:lagrange_multiplier) \n");
348 :
349 0 : 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 0 : m_spLagMultDisc->lagrange_multiplier(lagrangeMult, u, m_vActiveSetGlob, m_vActiveSubsets);
354 : }
355 0 : }
356 :
357 : template <typename TDomain, typename TAlgebra>
358 : template <typename TElem, typename TIterator>
359 : void ActiveSet<TDomain, TAlgebra>::lagrange_mat_inv_elem(TIterator iterBegin,
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>
432 : void ActiveSet<TDomain, TAlgebra>::lagrange_mat_inv(matrix_type& lagrangeMatInv)
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>
463 0 : void ActiveSet<TDomain, TAlgebra>::residual_lagrange_mult(vector_type& lagMult,
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 0 : if(m_vActiveSetGlob.size() != 0.0)
470 : {
471 0 : if (u.size() != lagMult.size())
472 0 : 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 0 : 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 0 : 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 0 : for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
494 0 : itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
495 : {
496 : // compute lagrange multiplier for active multiIndices
497 : // lagMult = rhs - Mat * u;
498 0 : 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 0 : }
514 :
515 :
516 : template <typename TDomain, typename TAlgebra>
517 : template <typename TElem, typename TIterator>
518 0 : bool ActiveSet<TDomain, TAlgebra>::check_conv_elem(TIterator iterBegin,
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 0 : = 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 0 : for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
541 : {
542 : TElem* elem = *iter;
543 :
544 : // reference object type
545 0 : 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 0 : int nCorner = (int)vCorner.size();
556 0 : for (int i = 0; i < nCorner; ++i)
557 0 : vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
558 :
559 0 : if (nCorner == dim)
560 0 : ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
561 : else
562 0 : 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 0 : locU.resize(ind); locObs.resize(indObs); locLambda.resize(ind);
569 :
570 : // read local values of u and lagrangeMult
571 0 : GetLocalVector(locU, u); GetLocalVector(locObs, *m_spObs);
572 0 : 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 0 : for(size_t fct = 0; fct < nrFctElem; ++fct)
580 : {
581 : size_t nrDoFsPerFctElem = ind.num_dof(fct);
582 0 : for(size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
583 : {
584 0 : for(int i = 0; i < dim; ++i)
585 0 : locUDof[i] = locU(i, dof);
586 0 : locUNormal = VecDot(locUDof, normal) / normOfNormal;
587 :
588 0 : gapValue = locUNormal - locObs(fct, dof);
589 0 : if (gapValue > 1e-06){
590 : // i.e.: m_spObs < u
591 : // constraint is violated
592 0 : return false;
593 : }
594 :
595 0 : kktcond = gapValue * locLambda(fct, dof);
596 0 : if (kktcond > 1e-06 || kktcond < -1e-06)
597 : return false;
598 :
599 : } // end(dof)
600 : } // end(fct)
601 : } // end(elem)
602 :
603 0 : return true;
604 0 : }
605 :
606 : template <typename TDomain, typename TAlgebra>
607 0 : bool ActiveSet<TDomain, TAlgebra>::check_conv(function_type& u, const function_type& lambda,
608 : const size_t step)
609 : {
610 : // ensure that at least one activeSet-iteration is performed
611 0 : 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 0 : if (m_vActiveSetGlob == m_vActiveSetGlobOld)
624 : {
625 : // check if constraint is fulfilled
626 : bool bConstraintViolated = false;
627 0 : for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
628 0 : activeSI != m_vActiveSubsets.end(); ++activeSI)
629 : {
630 0 : const int subsetDim = DimensionOfSubset(*m_spDD->subset_handler(), *activeSI);
631 0 : switch(subsetDim)
632 : {
633 : case 0:
634 : break;
635 0 : case 1:
636 0 : if (!check_conv_elem<RegularEdge>
637 0 : (m_spDD->template begin<RegularEdge>(*activeSI), m_spDD->template end<RegularEdge>(*activeSI),
638 : u, lambda))
639 : {bConstraintViolated = true;}
640 :
641 : break;
642 0 : case 2:
643 0 : if (!check_conv_elem<Triangle>
644 0 : (m_spDD->template begin<Triangle>(*activeSI),
645 : m_spDD->template end<Triangle>(*activeSI), u, lambda))
646 : {bConstraintViolated = true;}
647 :
648 0 : if (!check_conv_elem<Quadrilateral>
649 0 : (m_spDD->template begin<Quadrilateral>(*activeSI),
650 : m_spDD->template end<Quadrilateral>(*activeSI), u, lambda))
651 : {bConstraintViolated = true;}
652 :
653 : break;
654 0 : default:
655 0 : UG_THROW("ActiveSet::check_conv:"
656 : "SubsetDimension "<< subsetDim <<" (subset="<< *activeSI <<") not supported.");
657 : }
658 :
659 0 : 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>
673 0 : bool ActiveSet<TDomain, TAlgebra>::check_inequ(const matrix_type& mat,
674 : const vector_type& u,
675 : const vector_type& lambda,
676 : const vector_type& rhs)
677 : {
678 0 : if (u.size() != lambda.size())
679 0 : UG_THROW("Temporarily u and lambda need to be "
680 : "of same size in ActiveSet:check_inequ \n");
681 :
682 0 : SmartPtr<vector_type> spMat_u = u.clone_without_values();
683 0 : 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 0 : MatMult(*spMat_u, 1.0, mat, u);
692 : #endif
693 :
694 0 : for (size_t i = 0; i < u.size(); ++i)
695 : {
696 : //if (lambda[i] < -1e-06)
697 : // return false;
698 :
699 0 : (*spRes)[i] = (*spMat_u)[i] - rhs[i]; //+ lambda[i] - rhs[i]; //
700 0 : UG_LOG("lambda["<< i << "]: " << lambda[i] <<", res[" << i << "]: " << (*spRes)[i] << "\n");
701 : //if ((*spRes)[i] > 1e-06)
702 : // return false;
703 : }
704 :
705 0 : return true;
706 : }
707 :
708 : }; // namespace ug
709 :
710 : #endif /* ACTIVE_SET_IMPL_H_ */
|
