...
 
Commits (2)
......@@ -12,6 +12,30 @@ namespace Dune {
namespace Vem {
// template <class Matrix>
// struct ColumnVector
// {
// ColumnVector(Matrix &matrix, int col)
// : matrix_(matrix), col_(col) {}
// int size() const { return matrix_.rows(); }
// // typename Matrix::value_type& operator[](int i) {return matrix_[i][col_];}
// template <class Vector>
// ColumnVector &operator=(const Vector &v)
// {
// std::cout << "v size: " << v.size()
// << " size: " << size() << std::endl;
//
// assert( v.size() == size() );
// for ( std::size_t i = 0; i < size(); ++i)
// matrix_[i][col_] = v[i];
// }
// Matrix &matrix_;
// int col_;
// };
// template <class Matrix>
// ColumnVector<Matrix> columnVector(Matrix &matrix, int col)
// { return ColumnVector(matrix,col); }
template< class Matrix>
class LeastSquares {
public:
......@@ -33,7 +57,7 @@ namespace Dune {
template <class Vector>
Vector solve(const Vector &b, const Vector &d){
assert( b.size() == llsMatrix_.rows() );
Vector systemMultiply;
Vector systemMultiply, systemLagrange;
// std::cout << constraintMatrix_.size() << "is empty" << isEmpty(constraintMatrix_) << std::endl;
......@@ -50,11 +74,16 @@ namespace Dune {
Size systemMatrixDim = systemMatrixInv_.rows();
// since systemMatrix square, rows = cols
systemMultiply.resize(systemMatrixDim,0);
systemLagrange.resize(systemMatrixDim,0);
for (Size i = 0; i < systemMatrixDim; ++i)
for (Size j = 0; j < systemMatrixDim; ++j)
systemMultiply[i] += systemMatrixInv_[i][j] * systemVector[j];
systemLagrange[i] += systemMatrixInv_[i][j] * systemVector[j];
systemMultiply.resize(llsMatrix_.cols(),0);
// get rid of Lagrange multipliers
for( Size i = 0; i < systemMultiply.size(); ++i)
systemMultiply[i] = systemLagrange[i];
}
return systemMultiply;
}
......
......@@ -285,7 +285,7 @@ namespace Dune
ColumnVector &operator=(const Vector &v)
{
assert( v.size() == size() );
for (std::size_t i=0;i<size();++i)
for ( std::size_t i=0; i < size(); ++i)
matrix_[i][col_] = v[i];
}
Matrix &matrix_;
......@@ -295,6 +295,26 @@ namespace Dune
ColumnVector<Matrix> columnVector(Matrix &matrix, int col)
{ return ColumnVector(matrix,col); }
// L.solve(d,b,columnVector(valueProjection,beta));
template< class Matrix >
void printMatrix(const Matrix &A)
{
for(unsigned int i = 0; i < A.rows(); ++i) {
for (unsigned int j = 0; j < A.cols(); ++j) {
std::cout << A[i][j] << " ";
}
std::cout << std::endl;
}
}
template< class Vector>
void printVector(const Vector &x)
{
for(unsigned int i = 0; i < x.size(); ++i)
std::cout << x[i] << std::endl;
}
void buildProjections ();
// issue with making these const: use of delete default constructor in some python bindings...
......@@ -358,34 +378,39 @@ namespace Dune
size(std::max(polOrder-2,0));
#endif
std::cout << "size of spaces: "
<< numInnerShapeFunctions << " "
<< numShapeFunctions << " "
<< numGradShapeFunctions << " "
<< numHessShapeFunctions << std::endl;
// set up matrices used for constructing gradient, value, and edge projections
// Note: the code is set up with the assumption that the dofs suffice to compute the edge projection
DynamicMatrix< DomainFieldType > D, C, constraintValueProj, Hp, HpGrad, HpHess, HpInv, HpGradInv, HpHessInv;
DynamicMatrix< DomainType > R; // ,G //!!!
DynamicMatrix< typename Traits::ScalarBasisFunctionSetType::HessianMatrixType > P;
LeftPseudoInverse< DomainFieldType > pseudoInverse( numShapeFunctions );
// set up matrices used for constructing gradient, value, and edge projections
// Note: the code is set up with the assumption that the dofs suffice to compute the edge projection
DynamicMatrix< DomainFieldType > D, C, constraintValueProj, Hp, HpGrad, HpHess, HpInv, HpGradInv, HpHessInv;
DynamicMatrix< DomainType > R; // ,G //!!!
DynamicMatrix< typename Traits::ScalarBasisFunctionSetType::HessianMatrixType > P;
// these are the matrices we need to compute
valueProjections_.resize( agglomeration().size() );
jacobianProjections_.resize( agglomeration().size() );
hessianProjections_.resize( agglomeration().size() );
stabilizations_.resize( agglomeration().size() );
LeftPseudoInverse< DomainFieldType > pseudoInverse( numShapeFunctions );
// start iteration over all polygons
for( std::size_t agglomerate = 0; agglomerate < agglomeration().size(); ++agglomerate )
// these are the matrices we need to compute
valueProjections_.resize( agglomeration().size() );
jacobianProjections_.resize( agglomeration().size() );
hessianProjections_.resize( agglomeration().size() );
stabilizations_.resize( agglomeration().size() );
// start iteration over all polygons
for( std::size_t agglomerate = 0; agglomerate < agglomeration().size(); ++agglomerate )
{
const auto &bbox = blockMapper_.indexSet().boundingBox(agglomerate);
const std::size_t numDofs = blockMapper().numDofs( agglomerate );
const auto &bbox = blockMapper_.indexSet().boundingBox(agglomerate);
const std::size_t numDofs = blockMapper().numDofs( agglomerate );
std::cout << "Num dofs: "
<< numDofs << std::endl;
D.resize( numDofs, numShapeFunctions, 0 );
C.resize( numShapeFunctions, numDofs, 0 );
Hp.resize( numShapeFunctions, numShapeFunctions, 0 );
HpGrad.resize( numGradShapeFunctions, numGradShapeFunctions, 0 );
D.resize( numDofs, numShapeFunctions, 0 );
C.resize( numShapeFunctions, numDofs, 0 );
Hp.resize( numShapeFunctions, numShapeFunctions, 0 );
HpGrad.resize( numGradShapeFunctions, numGradShapeFunctions, 0 );
HpHess.resize( numHessShapeFunctions, numHessShapeFunctions, 0);
//!!! G.resize( numGradShapeFunctions, numGradShapeFunctions, DomainType(0) );
R.resize( numGradShapeFunctions, numDofs, DomainType(0) );
......@@ -448,15 +473,31 @@ namespace Dune
jacobianProjection[ alpha ].resize( numDofs, DomainType( 0 ) );
hessianProjection[ alpha ].resize( numDofs, 0 ); // typename hessianProjection[alpha]::value_type( 0 ) );
}
// valueProjection.resize( numDofs );
// for( std::size_t beta = 0; beta < numDofs; ++beta)
// valueProjection[ beta ].resize( numShapeFunctions, 0);
// type def for standard vector (to pick up re size for Hessian projection)
// need to resize Hessian projection
// re implementation of the value projection
auto leastSquaresMinimizer = LeastSquares( D, constraintValueProj );
DynamicVector< DomainFieldType > b( numDofs, 0 ), d( numInnerShapeFunctions, 0 );
std::vector< DomainFieldType > b( numDofs, 0 ), d( numInnerShapeFunctions, 0 );
// DynamicVector< DomainFieldType > solutionVector;
std::cout << "D: " << std::endl;
printMatrix(D);
std::cout << "constraints: " << std::endl;
printMatrix(constraintValueProj);
std::cout << "Hp: " << std::endl;
printMatrix(Hp);
std::cout << "I reached here" << std::endl;
for ( std::size_t beta = 0; beta < numDofs; ++beta ) {
auto colVecValueProjection = ColumnVector( valueProjection, beta );
auto colValueProjection = ColumnVector( valueProjection, beta );
// set up vectors b and d needed for least squares
b[ beta ] = 1;
......@@ -464,13 +505,32 @@ namespace Dune
d[ beta - numDofs + numInnerShapeFunctions] = 1;
}
colVecValueProjection = leastSquaresMinimizer.solve( b, d );
std::cout << "b: " << std::endl;
printVector(b);
std::cout << "d: " << std::endl;
printVector(d);
colValueProjection = leastSquaresMinimizer.solve( b, d );
// std::cout << "Solution vec " << std::endl;
// printVector(solutionVector);
// colValueProjection = solutionVector;
// re-set vectors b and d
d[ beta - numDofs + numInnerShapeFunctions] = 0;
b[ beta ] = 0;
}
std::cout << "Value projection " << std::endl;
for (std::size_t alpha=0; alpha<numShapeFunctions; ++alpha) {
for (std::size_t i=0; i<numDofs; ++i) {
std::cout << valueProjection[alpha][i] << " ";
}
std::cout << std::endl;
}
std::cout << "I reached here 2" << std::endl;
/////////////////////////////////////////
/////////////////////////////////////////
// !!! Original value projection implementation
......@@ -481,7 +541,7 @@ namespace Dune
if (1 || numInnerShapeFunctions > 0)
{
// modify C for inner dofs
std::size_t alpha=0;
// std::size_t alpha=0;
/////////////////////////////////////////
/////////////////////////////////////////
//// !!! Original value projection implementation
......
......@@ -106,12 +106,15 @@ int main( int argc, char **argv )
printVector(exactSoln);
auto leastSquaresMinimizer = Dune::Vem::LeastSquares(A,C);
Dune::DynamicVector< double > solutionVector;
for( unsigned int i = 0; i < valueProj.cols(); ++i){
auto colVecAdjustment = Dune::Vem::ColumnVector(valueProj,i);
auto columnValueProj = Dune::Vem::ColumnVector(valueProj,i);
// Dune::Vem::ColumnVector<Dune::DynamicMatrix<double>> columnVectorVP(valueProj,i);
leastSquaresMinimizer.solve(b,d,colVecAdjustment);
columnValueProj = leastSquaresMinimizer.solve(b,d);
// printVector(solutionVector);
// columnValueProj = (solutionVector);
}
std::cout << "value projection " << std::endl;
......@@ -217,7 +220,7 @@ int main( int argc, char **argv )
auto colVecAdjustment = Dune::Vem::ColumnVector(valueProj,i);
// Dune::Vem::ColumnVector<Dune::DynamicMatrix<double>> columnVectorVP(valueProj,i);
leastSquaresMinimizer.solve(b,d,colVecAdjustment);
leastSquaresMinimizer.solve(b,d);
}
std::cout << "value projection " << std::endl;
......@@ -330,7 +333,7 @@ int main( int argc, char **argv )
// Dune::Vem::ColumnVector<Dune::DynamicMatrix<double>> columnVectorVP(valueProj,i);
leastSquaresMinimizer.solve(b,d,colVecAdjustment);
leastSquaresMinimizer.solve(b,d);
}
std::cout << "value projection " << std::endl;
......