ugbase/lib_algebra/active_set/active_set_impl.h
/*
* Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
* Author: Raphael Prohl
*
* This file is part of UG4.
*
* UG4 is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License version 3 (as published by the
* Free Software Foundation) with the following additional attribution
* requirements (according to LGPL/GPL v3 §7):
*
* (1) The following notice must be displayed in the Appropriate Legal Notices
* of covered and combined works: "Based on UG4 (www.ug4.org/license)".
*
* (2) The following notice must be displayed at a prominent place in the
* terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
*
* (3) The following bibliography is recommended for citation and must be
* preserved in all covered files:
* "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
* parallel geometric multigrid solver on hierarchically distributed grids.
* Computing and visualization in science 16, 4 (2013), 151-164"
* "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
* flexible software system for simulating pde based models on high performance
* computers. Computing and visualization in science 16, 4 (2013), 165-179"
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*/
#ifndef ACTIVE_SET_IMPL_H_
#define ACTIVE_SET_IMPL_H_
#include "active_set.h"
#include "lib_disc/common/geometry_util.h"
#include "lib_disc/spatial_disc/disc_util/fe_geom.h"
using namespace std;
namespace ug {
template <int dim> struct face_type_traits
{
typedef void face_type0;
typedef void face_type1;
typedef void DimFEGeo;
};
template <> struct face_type_traits<1>
{
typedef ReferenceVertex face_type0;
typedef ReferenceVertex face_type1;
typedef DimFEGeometry<1, 1> DimFEGeo;
};
template <> struct face_type_traits<2>
{
typedef ReferenceEdge face_type0;
typedef ReferenceEdge face_type1;
typedef DimFEGeometry<2, 1> DimFEGeo;
};
template <> struct face_type_traits<3>
{
typedef ReferenceTriangle face_type0;
typedef ReferenceQuadrilateral face_type1;
typedef DimFEGeometry<3, 2> DimFEGeo;
};
template <typename TDomain, typename TAlgebra>
void ActiveSet<TDomain, TAlgebra>::prepare(function_type& u)
{
m_vActiveSetGlob.resize(0); m_vActiveSetGlobOld.resize(0);
m_spDD = u.dof_distribution();
m_spDom = u.domain();
}
/*template <typename TDomain, typename TAlgebra>
bool ActiveSet<TDomain, TAlgebra>::check_dist_to_obs(vector_type& u)
{
// STILL IN PROGRESS: u sollte hier reference-position + Startl�sung sein!
value_type dist;
bool geometry_cut_by_cons = false;
for(size_t i = 0; i < u.size(); i++)
{
UG_LOG("u(" << i << "):" << u[i] << "\n");
UG_LOG("m_spObs(" << i << "):" << (*m_spObs)[i] << "\n");
dist = (*m_spObs)[i] - u[i];
UG_LOG("dist:" << dist << "\n");
//TODO: anstatt u muss hier die geometrische Info einflie�en!
for (size_t fct = 0; fct < m_nrFcts; fct++)
{
if (BlockRef(dist,fct) < 0.0) // i.e.: u < m_spObs
{
geometry_cut_by_cons = true;
break;
}
}
if (geometry_cut_by_cons)
break;
}
return geometry_cut_by_cons;
}*/
// builds up a vector of global (dof,fct)-pairs, which are active (wrt. an obstacle-constraint)
template <typename TDomain, typename TAlgebra>
template <typename TElem, typename TIterator>
void ActiveSet<TDomain, TAlgebra>::active_index_elem(TIterator iterBegin,
TIterator iterEnd,
function_type& u,
function_type& rhs,
function_type& lagrangeMult)
{
// check if at least an element exists, else return
if(iterBegin == iterEnd) return;
int elemCounter = 0;
static const int dim = function_type::dim;
typedef typename function_type::domain_type domain_type;
typedef typename face_type_traits<dim>::face_type0 face_type0;
typedef typename face_type_traits<dim>::face_type1 face_type1;
// get position accessor
typename domain_type::position_accessor_type& aaPos
= u.domain()->position_accessor();
// storage for corner coordinates
vector<MathVector<dim> > vCorner;
vector<MathVector<dim> > vSideCoPos;
MathVector<dim> normal;
// local indices and local algebra
LocalIndices ind, indObs;
LocalVector locU, locLM, locObs;
number complementaryVal;
// TODO: it could be more efficient to store the active elements
// and its active local indices as well
// Loop over all elements on active subsets
for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
{
TElem* sideElem = *iter;
// reference object type
ReferenceObjectID roid = sideElem->reference_object_id();
const DimReferenceElement<dim-1>& rRefElem
= ReferenceElementProvider::get<dim-1>(roid);
// get corners of element
CollectCornerCoordinates(vCorner, *sideElem, aaPos);
// here the ordering of the corners in the reference element is exploited
// in order to compute the outer normal later on
int nCorner = (int)vCorner.size();
for (int i = 0; i < nCorner; ++i)
vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
if (nCorner == dim)
{
ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
UG_LOG("face_type0 \n");
}
else
{
ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
UG_LOG("face_type1 \n");
}
/*for (int i = 0; i < (int)vCorner.size(); ++i)
UG_LOG("sideCoPos: " << sideCoPos[i] << "\n");
UG_LOG("normal: " << normal << "\n");*/
elemCounter++;
// get global indices
u.indices(*iter, ind); (*m_spObs).indices(*iter, indObs);
// adapt local algebra
locU.resize(ind); locLM.resize(ind); locObs.resize(indObs);
// read local values of u and lagrangeMult
GetLocalVector(locU, u); GetLocalVector(locLM, lagrangeMult);
GetLocalVector(locObs, *m_spObs);
// loop over DoFs in element and store all activeMultiIndex-pairs in vector
size_t nrFctElem = ind.num_fct();
number normOfNormal = VecLength(normal);
MathVector<dim> locUDof;
for(size_t fct = 0; fct < nrFctElem; ++fct)
{
size_t nrDoFsPerFctElem = ind.num_dof(fct);
for(size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
{
for(int i = 0; i < dim; ++i)
locUDof[i] = locU(i, dof);
// locLM corresponds to the lagrange multiplier lambda,
// c.f. Hintermueller/Ito/Kunisch(2003)
number locUNormal = VecDot(locUDof, normal) / normOfNormal;
complementaryVal = locLM(fct ,dof) + locUNormal - locObs(fct, dof);
if (complementaryVal <= 1e-10){
locLM(fct,dof) = 0.0;
}
else
{
// create list of active global MultiIndex-pairs. Only those pairs should be attached
// which are not already a member of the 'activeSetGlob'-vector
bool bAlreadyActive = false;
for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
{
if ((*itActiveInd)[0] == ind.index(fct, dof)
&& (*itActiveInd)[1] == ind.comp(fct, dof))
bAlreadyActive = true;
}
if (!bAlreadyActive)
{
DoFIndex newActiveIndex(ind.index(fct, dof), ind.comp(fct, dof));
m_vActiveSetGlob.push_back(newActiveIndex);
BlockRef(rhs[ind.index(fct, dof)], ind.comp(fct, dof)) = locObs.value(fct,dof);
}
}
} // end(dof)
} // end(fct)
// send local to global lagrangeMult
//AddLocalVector(lagrangeMult, locLM);
UG_LOG("#activeDoFFctPairs global: " << m_vActiveSetGlob.size() << "\n");
} // end(elem)
UG_LOG("#elems: " << elemCounter << "\n");
}
// builds up a vector of global (dof,fct)-pairs, which are active (wrt. an obstacle-constraint)
// and sets dirichlet values in rhs for active indices
template <typename TDomain, typename TAlgebra>
bool ActiveSet<TDomain, TAlgebra>::active_index(function_type& u,
function_type& rhs, function_type& lagrangeMult, function_type& gap)
{
// note: the active-index search is restricted
// to constraint of the form u * n <= consGF
if(!m_bObs)
UG_THROW("No constraint set in ActiveSet \n");
if (u.num_indices() != rhs.num_indices() || u.num_indices() != lagrangeMult.num_indices() )
UG_THROW("GridFunctions u, rhs and lagrangeMult need to be defined "
"for the same domain and of the same size in 'ActiveSet:active_index' \n");
// remember old ActiveSet for convergence check
// TODO: avoid this vector copy; m_vActiveSetGlobOld really necessary?
m_vActiveSetGlobOld = m_vActiveSetGlob;
m_vActiveSetGlob.resize(0);
// 1.) get all subsets on which the 'gap'-gridfunction is defined!
// -> store them in vSubsetsOflagrangeMults
m_vActiveSubsets.resize(0);
//TODO: it is only necessary to loop over all BoundarySubsets!
for (int si = 0; si < m_spDD->num_subsets(); si++){
for (size_t fct = 0; fct < gap.num_fct(); fct++)
if (gap.is_def_in_subset(fct, si))
{
m_vActiveSubsets.push_back(si);
// 'break' is necessary to ensure that 'si' is
// only added once when several fcts of
// 'gap' are defined in subset 'si'!
break;
}
}
if (m_vActiveSubsets.size() == 0)
UG_LOG("No subsets chosen as possible active subsets. \n");
UG_LOG("#sizeOfvActiveSubsets: " << m_vActiveSubsets.size() << "\n");
// 2.) loop over all elements of the possible active subsets
for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
activeSI != m_vActiveSubsets.end(); ++activeSI)
{
//UG_LOG("activeSI: " << *activeSI << "\n");
const int subsetDim = DimensionOfSubset(*m_spDD->subset_handler(), *activeSI);
//UG_LOG("subsetDim: " << subsetDim << "\n");
// 3.) get localU out of u for each element and
// 4.) store the active global (dof,fct)-pair in m_vActiveSetGlob
switch(subsetDim)
{
case 0:
break;
case 1:
active_index_elem<RegularEdge>
(m_spDD->template begin<RegularEdge>(*activeSI),
m_spDD->template end<RegularEdge>(*activeSI), u, rhs, lagrangeMult);
break;
case 2:
active_index_elem<Triangle>
(m_spDD->template begin<Triangle>(*activeSI),
m_spDD->template end<Triangle>(*activeSI), u, rhs, lagrangeMult);
active_index_elem<Quadrilateral>
(m_spDD->template begin<Quadrilateral>(*activeSI),
m_spDD->template end<Quadrilateral>(*activeSI), u, rhs, lagrangeMult);
break;
default:
UG_THROW("ActiveSet::active_index:"
"SubsetDimension "<< subsetDim <<" (subset="<< *activeSI <<") not supported.");
}
}
if (m_vActiveSetGlob.size() > 0) return true;
else return false;
}
// sets a Dirichlet row for all active Indices
template <typename TDomain, typename TAlgebra>
void ActiveSet<TDomain, TAlgebra>::
set_dirichlet_rows(matrix_type& mat)
{
for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd){
SetDirichletRow(mat, *itActiveInd);
}
}
// computes the lagrange multiplier for a given discretization
template <typename TDomain, typename TAlgebra>
void ActiveSet<TDomain, TAlgebra>::lagrange_multiplier(function_type& lagrangeMult,
const function_type& u)
{
// check that lagrange multiplier disc is set
if (m_spLagMultDisc.invalid())
UG_THROW("No discretization set to compute the lagrange multiplier (in "
"ActiveSet:lagrange_multiplier) \n");
if(m_vActiveSetGlob.size() != 0.0)
{
UG_LOG("activeDoFs in ActiveSet:lagrange_multiplier " << m_vActiveSetGlob.size() << "\n");
//TODO: pass localActiveElemAndIndex here!!!
m_spLagMultDisc->lagrange_multiplier(lagrangeMult, u, m_vActiveSetGlob, m_vActiveSubsets);
}
}
template <typename TDomain, typename TAlgebra>
template <typename TElem, typename TIterator>
void ActiveSet<TDomain, TAlgebra>::lagrange_mat_inv_elem(TIterator iterBegin,
TIterator iterEnd, matrix_type& lagrangeMatInv)
{
typedef typename face_type_traits<dim>::face_type0 face_type0;
typedef typename face_type_traits<dim>::face_type1 face_type1;
typedef typename face_type_traits<dim>::DimFEGeo sideGeo;
typename TDomain::position_accessor_type& aaPos = m_spDom->position_accessor();
// some storage
MathVector<dim> normal;
vector<MathVector<dim> > vCorner;
vector<MathVector<dim> > vSideCoPos;
// local indices and local algebra
LocalIndices ind; LocalMatrix loclagrangeMatInv;
// Loop over all elements on active subsets
for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
{
TElem* elem = *iter;
// get global indices
m_spDD->indices(elem, ind);
loclagrangeMatInv.resize(ind, ind);
loclagrangeMatInv = 0.0;
// reference object type and geometry
ReferenceObjectID sideRoid = elem->reference_object_id();
sideGeo geo(sideRoid, 3, LFEID(LFEID::LAGRANGE, dim, 1));
// prepare geometry for type and order
try{
geo.update_local(sideRoid, LFEID(LFEID::LAGRANGE, dim, 1), 1);
}
UG_CATCH_THROW("ActiveSet::lagrange_mat_inv_elem:"
" Cannot update local values of finite element geometry.");
// get corners of element
CollectCornerCoordinates(vCorner, *elem, aaPos);
const DimReferenceElement<dim-1>& rRefElem
= ReferenceElementProvider::get<dim-1>(sideRoid);
int nCorner = (int)vCorner.size();
for (int i = 0; i < nCorner; ++i)
vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
if (nCorner == dim)
ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
else
ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
UG_LOG("normal in lagrange_mat_inv_elem: " << normal << "\n");
//number normOfNormal = VecLength(normal);
try{
geo.update(elem, &vSideCoPos[0], LFEID(LFEID::LAGRANGE, dim, 1), 1);
}
UG_CATCH_THROW("ActiveSet::lagrange_mat_inv_elem:"
" Cannot update finite element geometry.");
for (size_t ip = 0; ip < geo.num_ip(); ++ip)
for(size_t sh = 0; sh < geo.num_sh(); ++sh)
for (size_t i = 0; i < (size_t) dim; ++i)
loclagrangeMatInv(i, sh, i, sh) += 1.0/(geo.weight(ip) * geo.shape(ip, sh));
// * normal[i] / normOfNormal;
AddLocalMatrixToGlobal(lagrangeMatInv, loclagrangeMatInv);
}
}
template <typename TDomain, typename TAlgebra>
void ActiveSet<TDomain, TAlgebra>::lagrange_mat_inv(matrix_type& lagrangeMatInv)
{
for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
activeSI != m_vActiveSubsets.end(); ++activeSI)
{
UG_LOG("activeSI: " << *activeSI << "\n");
const int subsetDim = DimensionOfSubset(*m_spDD->subset_handler(), *activeSI);
UG_LOG("subsetDim: " << subsetDim << "\n");
switch(subsetDim)
{
case 0:
break;
case 1:
lagrange_mat_inv_elem<RegularEdge>
(m_spDD->template begin<RegularEdge>(*activeSI), m_spDD->template end<RegularEdge>(*activeSI), lagrangeMatInv);
break;
case 2:
lagrange_mat_inv_elem<Triangle>
(m_spDD->template begin<Triangle>(*activeSI), m_spDD->template end<Triangle>(*activeSI), lagrangeMatInv);
lagrange_mat_inv_elem<Quadrilateral>
(m_spDD->template begin<Quadrilateral>(*activeSI), m_spDD->template end<Quadrilateral>(*activeSI), lagrangeMatInv);
break;
default:
UG_THROW("ActiveSet::lagrange_mat_inv:"
"SubsetDimension "<<subsetDim<<" (subset="<<*activeSI<<") not supported.");
}
}
}
template <typename TDomain, typename TAlgebra>
void ActiveSet<TDomain, TAlgebra>::residual_lagrange_mult(vector_type& lagMult,
const matrix_type& mat,
const vector_type& u,
vector_type& rhs)
{
// only if some indices are active the lagrange multiplier is computed
if(m_vActiveSetGlob.size() != 0.0)
{
if (u.size() != lagMult.size())
UG_THROW("Temporarily u and lagMult need to be "
"of same size in ActiveSet:residual_lagrange_mult \n");
/*matrix_type lagrangeMatInv;
//SmartPtr<AssembledLinearOperator<algebra_type> > splagrangeMatInv;
ass_lagrangeMatInv(lagrangeMatInv);*/
SmartPtr<vector_type> spMat_u = u.clone_without_values();
(*spMat_u).resize(u.size());
#ifdef UG_PARALLEL
MatMultDirect(*spMat_u, 1.0, mat, u);
spMat_u->set_storage_type(u.get_storage_mask());
#else
MatMult(*spMat_u, 1.0, mat, u);
#endif
//SmartPtr<vector_type> spRes = u.clone_without_values();
//(*spRes).resize(u.size());
// loop MultiIndex-pairs in activeSet-vector
for (vector<DoFIndex>::iterator itActiveInd = m_vActiveSetGlob.begin();
itActiveInd < m_vActiveSetGlob.end(); ++itActiveInd)
{
// compute lagrange multiplier for active multiIndices
// lagMult = rhs - Mat * u;
DoFRef(lagMult, *itActiveInd) = DoFRef(rhs, *itActiveInd) - DoFRef((*spMat_u), *itActiveInd);
}
/*#ifdef UG_PARALLEL
MatMultDirect(lagMult, 1.0, lagrangeMatInv, *spRes);
#else
MatMult(lagMult, 1.0, lagrangeMatInv, *spRes);
#endif*/
UG_LOG("new lagMult-values computed \n");
//UG_LOG("rhs updated \n");
}
else{
UG_LOG("no active index in residual_lagrange_mult \n");
}
}
template <typename TDomain, typename TAlgebra>
template <typename TElem, typename TIterator>
bool ActiveSet<TDomain, TAlgebra>::check_conv_elem(TIterator iterBegin,
TIterator iterEnd, function_type& u, const function_type& lambda)
{
static const int dim = function_type::dim;
typedef typename function_type::domain_type domain_type;
typedef typename face_type_traits<dim>::face_type0 face_type0;
typedef typename face_type_traits<dim>::face_type1 face_type1;
// get position accessor
typename domain_type::position_accessor_type& aaPos
= u.domain()->position_accessor();
// storage for corner coordinates
vector<MathVector<dim> > vCorner;
vector<MathVector<dim> > vSideCoPos;
MathVector<dim> normal;
// local indices and local algebra
LocalIndices ind, indObs;
LocalVector locU, locObs, locLambda;
// Loop over all elements on active subsets
for(TIterator iter = iterBegin; iter != iterEnd; ++iter)
{
TElem* elem = *iter;
// reference object type
ReferenceObjectID roid = elem->reference_object_id();
const DimReferenceElement<dim-1>& rRefElem
= ReferenceElementProvider::get<dim-1>(roid);
// get corners of element
CollectCornerCoordinates(vCorner, *elem, aaPos);
// here the ordering of the corners in the reference element is exploited
// in order to compute the outer normal later on
int nCorner = (int)vCorner.size();
for (int i = 0; i < nCorner; ++i)
vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, 0, 0, i)]);
if (nCorner == dim)
ElementNormal<face_type0, dim>(normal, &vSideCoPos[0]);
else
ElementNormal<face_type1, dim>(normal, &vSideCoPos[0]);
// get global indices
u.indices(*iter, ind); (*m_spObs).indices(*iter, indObs);
// adapt local algebra
locU.resize(ind); locObs.resize(indObs); locLambda.resize(ind);
// read local values of u and lagrangeMult
GetLocalVector(locU, u); GetLocalVector(locObs, *m_spObs);
GetLocalVector(locLambda, lambda);
size_t nrFctElem = ind.num_fct();
number gapValue, locUNormal, kktcond;
number normOfNormal = VecLength(normal);
MathVector<dim> locUDof;
for(size_t fct = 0; fct < nrFctElem; ++fct)
{
size_t nrDoFsPerFctElem = ind.num_dof(fct);
for(size_t dof = 0; dof < nrDoFsPerFctElem; ++dof)
{
for(int i = 0; i < dim; ++i)
locUDof[i] = locU(i, dof);
locUNormal = VecDot(locUDof, normal) / normOfNormal;
gapValue = locUNormal - locObs(fct, dof);
if (gapValue > 1e-06){
// i.e.: m_spObs < u
// constraint is violated
return false;
}
kktcond = gapValue * locLambda(fct, dof);
if (kktcond > 1e-06 || kktcond < -1e-06)
return false;
} // end(dof)
} // end(fct)
} // end(elem)
return true;
}
template <typename TDomain, typename TAlgebra>
bool ActiveSet<TDomain, TAlgebra>::check_conv(function_type& u, const function_type& lambda,
const size_t step)
{
// ensure that at least one activeSet-iteration is performed
if (step <= 1)
return false;
// NOW TWO CHECKS WILL BE PERFORMED TO ENSURE CONVERGENCE:
// 1. Did some multiIndices change from 'active' to 'inactive' or vice versa
// in the last iteration-step?
// 2. Is the constraint violated for any DoFIndex?
UG_LOG(m_vActiveSetGlob.size() << " indices are active at the begin "
"of step " << step << " ! \n");
// check if activeSet has changed
if (m_vActiveSetGlob == m_vActiveSetGlobOld)
{
// check if constraint is fulfilled
bool bConstraintViolated = false;
for (vector<int>::iterator activeSI = m_vActiveSubsets.begin();
activeSI != m_vActiveSubsets.end(); ++activeSI)
{
const int subsetDim = DimensionOfSubset(*m_spDD->subset_handler(), *activeSI);
switch(subsetDim)
{
case 0:
break;
case 1:
if (!check_conv_elem<RegularEdge>
(m_spDD->template begin<RegularEdge>(*activeSI), m_spDD->template end<RegularEdge>(*activeSI),
u, lambda))
{bConstraintViolated = true;}
break;
case 2:
if (!check_conv_elem<Triangle>
(m_spDD->template begin<Triangle>(*activeSI),
m_spDD->template end<Triangle>(*activeSI), u, lambda))
{bConstraintViolated = true;}
if (!check_conv_elem<Quadrilateral>
(m_spDD->template begin<Quadrilateral>(*activeSI),
m_spDD->template end<Quadrilateral>(*activeSI), u, lambda))
{bConstraintViolated = true;}
break;
default:
UG_THROW("ActiveSet::check_conv:"
"SubsetDimension "<< subsetDim <<" (subset="<< *activeSI <<") not supported.");
}
if (bConstraintViolated)
return false;
}
// activeSet remains unchanged & constraint is fulfilled for all indices
return true;
}
else{
return false;
}
}
template <typename TDomain, typename TAlgebra>
bool ActiveSet<TDomain, TAlgebra>::check_inequ(const matrix_type& mat,
const vector_type& u,
const vector_type& lambda,
const vector_type& rhs)
{
if (u.size() != lambda.size())
UG_THROW("Temporarily u and lambda need to be "
"of same size in ActiveSet:check_inequ \n");
SmartPtr<vector_type> spMat_u = u.clone_without_values();
SmartPtr<vector_type> spRes = u.clone_without_values();
(*spMat_u).resize(u.size()); (*spRes).resize(u.size());
#ifdef UG_PARALLEL
MatMultDirect(*spMat_u, 1.0, mat, u);
spMat_u->set_storage_type(u.get_storage_mask());
spRes->set_storage_type(u.get_storage_mask());
#else
MatMult(*spMat_u, 1.0, mat, u);
#endif
for (size_t i = 0; i < u.size(); ++i)
{
//if (lambda[i] < -1e-06)
// return false;
(*spRes)[i] = (*spMat_u)[i] - rhs[i]; //+ lambda[i] - rhs[i]; //
UG_LOG("lambda["<< i << "]: " << lambda[i] <<", res[" << i << "]: " << (*spRes)[i] << "\n");
//if ((*spRes)[i] > 1e-06)
// return false;
}
return true;
}
}; // namespace ug
#endif /* ACTIVE_SET_IMPL_H_ */