33 #ifndef ACTIVE_SET_IMPL_H_
34 #define ACTIVE_SET_IMPL_H_
72 template <
typename TDomain,
typename TAlgebra>
75 m_vActiveSetGlob.resize(0); m_vActiveSetGlobOld.resize(0);
113 template <
typename TDomain,
typename TAlgebra>
114 template <
typename TElem,
typename TIterator>
122 if(iterBegin == iterEnd)
return;
132 typename domain_type::position_accessor_type& aaPos
133 = u.
domain()->position_accessor();
136 vector<MathVector<dim> > vCorner;
137 vector<MathVector<dim> > vSideCoPos;
150 for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
152 TElem* sideElem = *iter;
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)]);
171 ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
176 ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
187 u.
indices(*iter, ind); (*m_spObs).indices(*iter, indObs);
197 size_t nrFctElem = ind.
num_fct();
201 for(
size_t fct = 0; fct < nrFctElem; ++fct)
203 size_t nrDoFsPerFctElem = ind.
num_dof(fct);
204 for(
size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
206 for(
int i = 0; i <
dim; ++i)
207 locUDof[i] = locU(i, dof);
211 number locUNormal =
VecDot(locUDof, normal) / normOfNormal;
212 complementaryVal = locLM(fct ,dof) + locUNormal - locObs(fct, dof);
214 if (complementaryVal <= 1e-10){
215 locLM(fct,dof) = 0.0;
221 bool bAlreadyActive =
false;
222 for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
223 itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
225 if ((*itActiveInd)[0] == ind.
index(fct, dof)
226 && (*itActiveInd)[1] == ind.
comp(fct, dof))
227 bAlreadyActive =
true;
233 m_vActiveSetGlob.push_back(newActiveIndex);
244 UG_LOG(
"#activeDoFFctPairs global: " << m_vActiveSetGlob.size() <<
"\n");
246 UG_LOG(
"#elems: " << elemCounter <<
"\n");
252 template <
typename TDomain,
typename TAlgebra>
259 UG_THROW(
"No constraint set in ActiveSet \n");
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");
267 m_vActiveSetGlobOld = m_vActiveSetGlob;
268 m_vActiveSetGlob.resize(0);
272 m_vActiveSubsets.resize(0);
274 for (
int si = 0; si < m_spDD->num_subsets(); si++){
275 for (
size_t fct = 0; fct < gap.
num_fct(); fct++)
278 m_vActiveSubsets.push_back(si);
286 if (m_vActiveSubsets.size() == 0)
287 UG_LOG(
"No subsets chosen as possible active subsets. \n");
289 UG_LOG(
"#sizeOfvActiveSubsets: " << m_vActiveSubsets.size() <<
"\n");
292 for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
293 activeSI != m_vActiveSubsets.end(); ++activeSI)
306 active_index_elem<RegularEdge>
307 (m_spDD->template begin<RegularEdge>(*activeSI),
308 m_spDD->template end<RegularEdge>(*activeSI), u, rhs, lagrangeMult);
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);
320 "SubsetDimension "<< subsetDim <<
" (subset="<< *activeSI <<
") not supported.");
324 if (m_vActiveSetGlob.size() > 0)
return true;
329 template <
typename TDomain,
typename TAlgebra>
333 for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
334 itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd){
340 template <
typename TDomain,
typename TAlgebra>
345 if (m_spLagMultDisc.invalid())
346 UG_THROW(
"No discretization set to compute the lagrange multiplier (in "
347 "ActiveSet:lagrange_multiplier) \n");
349 if(m_vActiveSetGlob.size() != 0.0)
351 UG_LOG(
"activeDoFs in ActiveSet:lagrange_multiplier " << m_vActiveSetGlob.size() <<
"\n");
353 m_spLagMultDisc->lagrange_multiplier(lagrangeMult, u, m_vActiveSetGlob, m_vActiveSubsets);
357 template <
typename TDomain,
typename TAlgebra>
358 template <
typename TElem,
typename TIterator>
365 typename TDomain::position_accessor_type& aaPos = m_spDom->position_accessor();
369 vector<MathVector<dim> > vCorner;
370 vector<MathVector<dim> > vSideCoPos;
376 for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
381 m_spDD->indices(elem, ind);
383 loclagrangeMatInv.
resize(ind, ind);
384 loclagrangeMatInv = 0.0;
388 sideGeo geo(sideRoid, 3,
LFEID(LFEID::LAGRANGE,
dim, 1));
392 geo.update_local(sideRoid,
LFEID(LFEID::LAGRANGE,
dim, 1), 1);
395 " Cannot update local values of finite element geometry.");
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)]);
408 ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
410 ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
412 UG_LOG(
"normal in lagrange_mat_inv_elem: " << normal <<
"\n");
416 geo.update(elem, &vSideCoPos[0],
LFEID(LFEID::LAGRANGE,
dim, 1), 1);
419 " Cannot update finite element geometry.");
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));
431 template <
typename TDomain,
typename TAlgebra>
434 for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
435 activeSI != m_vActiveSubsets.end(); ++activeSI)
437 UG_LOG(
"activeSI: " << *activeSI <<
"\n");
439 UG_LOG(
"subsetDim: " << subsetDim <<
"\n");
446 lagrange_mat_inv_elem<RegularEdge>
447 (m_spDD->template begin<RegularEdge>(*activeSI), m_spDD->template end<RegularEdge>(*activeSI), lagrangeMatInv);
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);
456 UG_THROW(
"ActiveSet::lagrange_mat_inv:"
457 "SubsetDimension "<<subsetDim<<
" (subset="<<*activeSI<<
") not supported.");
462 template <
typename TDomain,
typename TAlgebra>
469 if(m_vActiveSetGlob.size() != 0.0)
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");
480 (*spMat_u).resize(u.size());
484 spMat_u->set_storage_type(u.get_storage_mask());
486 MatMult(*spMat_u, 1.0, mat, u);
493 for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
494 itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
498 DoFRef(lagMult, *itActiveInd) =
DoFRef(rhs, *itActiveInd) -
DoFRef((*spMat_u), *itActiveInd);
507 UG_LOG(
"new lagMult-values computed \n");
511 UG_LOG(
"no active index in residual_lagrange_mult \n");
516 template <
typename TDomain,
typename TAlgebra>
517 template <
typename TElem,
typename TIterator>
527 typename domain_type::position_accessor_type& aaPos
528 = u.
domain()->position_accessor();
531 vector<MathVector<dim> > vCorner;
532 vector<MathVector<dim> > vSideCoPos;
540 for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
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)]);
560 ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
562 ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
565 u.
indices(*iter, ind); (*m_spObs).indices(*iter, indObs);
574 size_t nrFctElem = ind.
num_fct();
575 number gapValue, locUNormal, kktcond;
579 for(
size_t fct = 0; fct < nrFctElem; ++fct)
581 size_t nrDoFsPerFctElem = ind.
num_dof(fct);
582 for(
size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
584 for(
int i = 0; i <
dim; ++i)
585 locUDof[i] = locU(i, dof);
586 locUNormal =
VecDot(locUDof, normal) / normOfNormal;
588 gapValue = locUNormal - locObs(fct, dof);
589 if (gapValue > 1e-06){
595 kktcond = gapValue * locLambda(fct, dof);
596 if (kktcond > 1e-06 || kktcond < -1e-06)
606 template <
typename TDomain,
typename TAlgebra>
619 UG_LOG(m_vActiveSetGlob.size() <<
" indices are active at the begin "
620 "of step " << step <<
" ! \n");
623 if (m_vActiveSetGlob == m_vActiveSetGlobOld)
626 bool bConstraintViolated =
false;
627 for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
628 activeSI != m_vActiveSubsets.end(); ++activeSI)
636 if (!check_conv_elem<RegularEdge>
637 (m_spDD->template begin<RegularEdge>(*activeSI), m_spDD->template end<RegularEdge>(*activeSI),
639 {bConstraintViolated =
true;}
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;}
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;}
656 "SubsetDimension "<< subsetDim <<
" (subset="<< *activeSI <<
") not supported.");
659 if (bConstraintViolated)
672 template <
typename TDomain,
typename TAlgebra>
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");
684 (*spMat_u).resize(u.size()); (*spRes).resize(u.size());
688 spMat_u->set_storage_type(u.get_storage_mask());
689 spRes->set_storage_type(u.get_storage_mask());
691 MatMult(*spMat_u, 1.0, mat, u);
694 for (
size_t i = 0; i < u.size(); ++i)
699 (*spRes)[i] = (*spMat_u)[i] - rhs[i];
700 UG_LOG(
"lambda["<< i <<
"]: " << lambda[i] <<
", res[" << i <<
"]: " << (*spRes)[i] <<
"\n");
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
#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
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