Commit 22e128c0 authored by Alice Hodson's avatar Alice Hodson

starting to integrate lls into the value projection

parent 17e44512
......@@ -8,8 +8,6 @@
#include <assert.h>
//using namespace std;
namespace Dune {
namespace Vem {
......@@ -37,7 +35,11 @@ namespace Dune {
// template < class Field >
LeastSquares(const Matrix &llsMatrix, const Matrix &constraintMatrix)
: llsMatrix_(llsMatrix), constraintMatrix_(constraintMatrix),
systemMatrixInv_(matrixSetUp())
systemMatrixInv_(matrixSetUp(llsMatrix_, constraintMatrix_))
{
}
LeastSquares(const Matrix &llsMatrix)
: llsMatrix_(llsMatrix), constraintMatrix_(emptyMatrix()), systemMatrixInv_(matrixSetUp(llsMatrix_))
{
}
LeastSquares(const LeastSquares &source) = delete;
......@@ -49,11 +51,14 @@ namespace Dune {
Vector systemMultiply;
if ( constraintMatrix_.size() == 0 ){
std::cout << constraintMatrix_.size() << "is empty" << isEmpty(constraintMatrix_) << std::endl;
if ( isEmpty(constraintMatrix_) ){
systemMultiply.resize(llsMatrix_.cols());
systemMatrixInv_.mv(b,systemMultiply);
}
else {
std::cout << "constraint rows" << constraintMatrix_.rows() << std::endl;
assert(d.size() == constraintMatrix_.rows());
Vector systemVector = vectorSetUp(b, d);
......@@ -111,28 +116,49 @@ namespace Dune {
private:
const Matrix &llsMatrix_;
const Matrix &constraintMatrix_;
const Matrix constraintMatrix_;
const Matrix systemMatrixInv_;
Matrix emptyMatrix()
{
// return matrix with no size
Matrix A;
std::cout << "size of empty matrix " << A.size() << std::endl;
return A;
}
bool isEmpty(const Matrix &A)
{
if ( A.size() == 0 ) {
return true;
}
return false;
}
// template <class Field>
Matrix matrixSetUp()
Matrix matrixSetUp(const Matrix &llsMatrix_)
{
if ( constraintMatrix_.size() == 0) {
// this needs changing
LeftPseudoInverse< double > pseudoInverse( llsMatrix_.cols() );
// !!! this needs changing
LeftPseudoInverse< double > 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;
// 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;
llsMatrixPseudoInv.resize( llsMatrix_.cols(), llsMatrix_.rows() );
pseudoInverse( llsMatrix_, llsMatrixPseudoInv);
pseudoInverse( llsMatrix_, llsMatrixPseudoInv);
std::cout << "pseudo Inv of A " << std::endl;
printMatrix(llsMatrixPseudoInv);
std::cout << "pseudo Inv of A " << std::endl;
printMatrix(llsMatrixPseudoInv);
return llsMatrixPseudoInv;
}
return llsMatrixPseudoInv;
Matrix matrixSetUp(const Matrix &llsMatrix_, const Matrix &constraintMatrix_)
{
if ( isEmpty(constraintMatrix_) ) {
return matrixSetUp(llsMatrix_);
}
else {
// construct the matrix [2A^T*A C^T ; C 0] needed for least squares solution
......@@ -210,9 +236,13 @@ namespace Dune {
};
template <class Matrix>
template <class Matrix>
LeastSquares<Matrix> leastSquares(const Matrix &llsMatrix, const Matrix &constraintMatrix)
{ return LeastSquares<Matrix>(llsMatrix,constraintMatrix); }
template <class Matrix>
LeastSquares<Matrix> leastSquares(const Matrix &llsMatrix)
{ return LeastSquares<Matrix>(llsMatrix); }
} // namespace Vem
} // namespace Dune
......
......@@ -22,6 +22,7 @@
#include <dune/vem/agglomeration/dofmapper.hh>
#include <dune/vem/misc/compatibility.hh>
#include <dune/vem/misc/pseudoinverse.hh>
#include <dune/vem/misc/leastSquares.hh>
#include <dune/vem/agglomeration/dgspace.hh>
#include <dune/vem/space/basisfunctionset.hh>
#include <dune/vem/space/interpolation.hh>
......@@ -356,13 +357,13 @@ namespace Dune
// 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, Hp, HpGrad, HpHess, HpInv, HpGradInv, HpHessInv;
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 );
// these are the matrices we need to comput e
// these are the matrices we need to compute
valueProjections_.resize( agglomeration().size() );
jacobianProjections_.resize( agglomeration().size() );
hessianProjections_.resize( agglomeration().size() );
......@@ -382,6 +383,7 @@ namespace Dune
//!!! G.resize( numGradShapeFunctions, numGradShapeFunctions, DomainType(0) );
R.resize( numGradShapeFunctions, numDofs, DomainType(0) );
P.resize( numHessShapeFunctions, numDofs, 0);
constraintValueProj.resize( numInnerShapeFunctions, numShapeFunctions, 0 );
// iterate over the triangles of this polygon
for( const ElementSeedType &entitySeed : entitySeeds[ agglomerate ] )
......@@ -408,7 +410,9 @@ namespace Dune
HpGrad[alpha][beta] += phi[0]*psi[0]*weight;
if (alpha<numHessShapeFunctions && beta <numHessShapeFunctions)
HpHess[alpha][beta] += phi[0]*psi[0]*weight;
} );
if ( alpha<numInnerShapeFunctions )
constraintValueProj[alpha][beta] += phi[0]*psi[0]*weight;
} );
#if 0 // !!!!
if (alpha<numGradShapeFunctions)
shapeFunctionSet.jacobianEach( quadrature[ qp ], [ & ] ( std::size_t beta, typename ScalarShapeFunctionSetType::JacobianRangeType psi ) {
......@@ -438,18 +442,51 @@ namespace Dune
// type def for standard vector (to pick up re size for Hessian projection)
// need to resize Hessian projection
pseudoInverse( D, valueProjection );
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 ) {
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;
}
leastSquaresMinimizer.solve( b, d, colVecValueProjection );
// re-set vectors b and d
d[numDofs - beta ] = 0;
b[ beta ] = 0;
}
/////////////////////////////////////////
/////////////////////////////////////////
// !!! Original value projection implementation
// pseudoInverse( D, valueProjection );
/////////////////////////////////////////
/////////////////////////////////////////
if (1 || numInnerShapeFunctions > 0)
{
// modify C for inner dofs
std::size_t alpha=0;
for (; alpha<numInnerShapeFunctions; ++alpha)
C[alpha][alpha+numDofs-numInnerShapeFunctions] = H0;
for (; alpha<numShapeFunctions; ++alpha)
for (std::size_t i=0; i<numDofs; ++i)
for (std::size_t beta=0; beta<numShapeFunctions; ++beta)
C[alpha][i] += Hp[alpha][beta]*valueProjection[beta][i];
/////////////////////////////////////////
/////////////////////////////////////////
// !!! Original value projection implementation
// for (; alpha<numInnerShapeFunctions; ++alpha)
// C[alpha][alpha+numDofs-numInnerShapeFunctions] = H0;
// for (; alpha<numShapeFunctions; ++alpha)
// for (std::size_t i=0; i<numDofs; ++i)
// for (std::size_t beta=0; beta<numShapeFunctions; ++beta)
// C[alpha][i] += Hp[alpha][beta]*valueProjection[beta][i];
/////////////////////////////////////////
/////////////////////////////////////////
HpInv = Hp;
HpInv.invert();
......@@ -516,20 +553,25 @@ namespace Dune
// for (std::size_t beta=0; beta<numInnerShapeFunctions; ++beta)
// R[alpha][beta+numDofs-numInnerShapeFunctions].axpy(-H0, G[beta][alpha]);
/////////////////////////////////////////
/////////////////////////////////////////
// !!! Original value projection implementation
// now compute projection by multiplying with inverse mass matrix
for (std::size_t alpha=0; alpha<numShapeFunctions; ++alpha)
{
for (std::size_t i=0; i<numDofs; ++i)
{
valueProjection[alpha][i] = 0;
for (std::size_t beta=0; beta<numShapeFunctions; ++beta)
valueProjection[alpha][i] += HpInv[alpha][beta]*C[beta][i];
//!!!!
// if (alpha<numGradShapeFunctions)
// for (std::size_t beta=0; beta<numGradShapeFunctions; ++beta)
// jacobianProjection[alpha][i].axpy(HpGradInv[alpha][beta],R[beta][i]);
}
}
// for (std::size_t alpha=0; alpha<numShapeFunctions; ++alpha)
// {
// for (std::size_t i=0; i<numDofs; ++i)
// {
// valueProjection[alpha][i] = 0;
// for (std::size_t beta=0; beta<numShapeFunctions; ++beta)
// valueProjection[alpha][i] += HpInv[alpha][beta]*C[beta][i];
// //!!!!
// // if (alpha<numGradShapeFunctions)
// // for (std::size_t beta=0; beta<numGradShapeFunctions; ++beta)
// // jacobianProjection[alpha][i].axpy(HpGradInv[alpha][beta],R[beta][i]);
// }
// }
/////////////////////////////////////////
/////////////////////////////////////////
} // have some inner moments
//////////////////////////////////////////////////////////////////////////
......
......@@ -3,6 +3,16 @@
#include <dune/common/dynvector.hh>
#include <dune/common/dynmatrix.hh>
template <class Matrix>
Matrix emptyMatrix()
{
// return matrix with no size
Matrix A;
// std::cout << "size of empty matrix " << A.size() << std::endl;
return A;
}
template< class Matrix >
void printMatrix(const Matrix &A)
{
......@@ -144,6 +154,9 @@ int main( int argc, char **argv )
Dune::DynamicMatrix< double > valueProj;
Dune::DynamicVector< double > exactSoln(3,0);
std::cout << "Sie of C " << C.size() << std::endl;
const int matrixDim = 3;
valueProj.resize( 3, 2, 0);
......@@ -196,7 +209,7 @@ int main( int argc, char **argv )
std::cout << "Vector exact solution: " << std::endl;
printVector(exactSoln);
auto leastSquaresMinimizer = Dune::Vem::LeastSquares(A,C);
auto leastSquaresMinimizer = Dune::Vem::LeastSquares(A);
std::cout << "I reached here" << std::endl;
......
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