Commit 8e7e5c10 authored by Andreas Dedner's avatar Andreas Dedner

fixed some bugs but still not working

parent 22e128c0
......@@ -12,25 +12,11 @@ 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_];}
Matrix &matrix_;
int col_;
};
template <class Matrix>
ColumnVector<Matrix> columnVector(Matrix &matrix, int col)
{ return ColumnVector(matrix,col); }
template< class Matrix>
class LeastSquares {
public:
typedef typename Matrix::size_type Size;
typedef typename Matrix::value_type Field;
// template < class Field >
LeastSquares(const Matrix &llsMatrix, const Matrix &constraintMatrix)
......@@ -44,11 +30,9 @@ namespace Dune {
}
LeastSquares(const LeastSquares &source) = delete;
template <class Vector, class Solution>
void solve(const Vector &b, const Vector &d, Solution &solution){
template <class Vector>
Vector solve(const Vector &b, const Vector &d){
assert( b.size() == llsMatrix_.rows() );
assert( solution.size() == llsMatrix_.cols() );
Vector systemMultiply;
std::cout << constraintMatrix_.size() << "is empty" << isEmpty(constraintMatrix_) << std::endl;
......@@ -72,15 +56,7 @@ namespace Dune {
for (Size j = 0; j < systemMatrixDim; ++j)
systemMultiply[i] += systemMatrixInv_[i][j] * systemVector[j];
}
for (Size i = 0; i < solution.size(); ++i)
solution[i] = systemMultiply[i];
}
template <class Vector>
Vector solve(const Vector &b, const Vector &d ){
Vector solution(llsMatrix_.cols());
solve(b,d,solution);
return solution;
return systemMultiply;
}
template <class Vector>
......@@ -97,6 +73,7 @@ namespace Dune {
void printMatrix(const Matrix &A)
{
return;
for(unsigned int i = 0; i < A.rows(); ++i) {
for (unsigned int j = 0; j < A.cols(); ++j) {
std::cout << A[i][j];
......@@ -108,6 +85,7 @@ namespace Dune {
template< class Vector>
void printVector(const Vector &x)
{
return;
for(unsigned int i = 0; i < x.size(); ++i)
std::cout << x[i] << std::endl;
}
......@@ -129,24 +107,20 @@ namespace Dune {
bool isEmpty(const Matrix &A)
{
if ( A.size() == 0 ) {
return true;
}
return false;
return ( A.size() == 0 );
}
// template <class Field>
Matrix matrixSetUp(const Matrix &llsMatrix_)
{
// !!! this needs changing
LeftPseudoInverse< double > pseudoInverse( llsMatrix_.cols() );
LeftPseudoInverse< Field > pseudoInverse( llsMatrix_.cols() );
// LeftPseudoInverse< Field > pseudoInverse( llsMatrix_.cols() );
// no constraints in this case and so form pseudo inverse
std::cout << "Matrix C has no size" << std::endl;
Matrix llsMatrixPseudoInv;
llsMatrixPseudoInv.resize( llsMatrix_.cols(), llsMatrix_.rows() );
Matrix llsMatrixPseudoInv( llsMatrix_.cols(), llsMatrix_.rows() );
pseudoInverse( llsMatrix_, llsMatrixPseudoInv);
......
......@@ -279,8 +279,15 @@ namespace Dune
{
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_];}
int size() const { return matrix_.size(); }
// typename Matrix::value_type& operator[](int i) {return matrix_[i][col_];}
template <class Vector>
ColumnVector &operator=(const Vector &v)
{
assert( v.size() == size() );
for (std::size_t i=0;i<size();++i)
matrix_[i][col_] = v[i];
}
Matrix &matrix_;
int col_;
};
......@@ -433,35 +440,37 @@ namespace Dune
auto &valueProjection = valueProjections_[ agglomerate ];
auto &jacobianProjection = jacobianProjections_[ agglomerate ];
auto &hessianProjection = hessianProjections_[ agglomerate ];
valueProjection.resize( numShapeFunctions );
jacobianProjection.resize( numShapeFunctions );
hessianProjection.resize( numShapeFunctions );
for( std::size_t alpha = 0; alpha < numShapeFunctions; ++alpha ){
valueProjection[ alpha ].resize( numDofs, 0 );
jacobianProjection[ alpha ].resize( numDofs, DomainType( 0 ) );
hessianProjection[ alpha ].resize( numDofs, 0 ); // typename hessianProjection[alpha]::value_type( 0 ) );
}
// type def for standard vector (to pick up re size for Hessian projection)
// need to resize Hessian projection
printMatrix(D);
printMatrix(constraintValueProj);
// printMatrix(D);
// printMatrix(constraintValueProj);
// re implementation of the value projection
auto leastSquaresMinimizer = LeastSquares( D, constraintValueProj );
DynamicVector< DomainFieldType > b( numDofs, 0 ), d( numInnerShapeFunctions, 0 );
for ( std::size_t beta = 0; beta < numShapeFunctions; ++beta ) {
for ( std::size_t beta = 0; beta < numDofs; ++beta ) {
auto colVecValueProjection = ColumnVector( valueProjection, beta );
// set up vectors b and d needed for least squares
b[ beta ] = 1;
if( beta > numDofs - numInnerShapeFunctions ){
d[ numDofs - beta ] = 1;
if( beta >= numDofs - numInnerShapeFunctions ){
d[ beta - numDofs + numInnerShapeFunctions] = 1;
}
leastSquaresMinimizer.solve( b, d, colVecValueProjection );
colVecValueProjection = leastSquaresMinimizer.solve( b, d );
// re-set vectors b and d
d[numDofs - beta ] = 0;
d[ beta - numDofs + numInnerShapeFunctions] = 0;
b[ beta ] = 0;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment