Commit 460bf883 authored by Andreas Dedner's avatar Andreas Dedner

fixed a bug in the lls constraint computation of the value projection

parent e891b47a
......@@ -49,7 +49,7 @@ namespace Dune {
{
}
LeastSquares(const Matrix &llsMatrix)
: llsMatrix_(llsMatrix), constraintMatrix_(emptyMatrix()), systemMatrixInv_(matrixSetUp(llsMatrix_))
: llsMatrix_(llsMatrix), constraintMatrix_(), systemMatrixInv_(matrixSetUp(llsMatrix_))
{
}
LeastSquares(const LeastSquares &source) = delete;
......@@ -82,6 +82,7 @@ namespace Dune {
systemMultiply.resize(llsMatrix_.cols(),0);
// get rid of Lagrange multipliers
// TODO? avoid copy by cutting of in operator= of ColVec
for( Size i = 0; i < systemMultiply.size(); ++i)
systemMultiply[i] = systemLagrange[i];
}
......@@ -123,17 +124,10 @@ namespace Dune {
private:
const Matrix &llsMatrix_;
// TODO: avoid copy of constraintMatrix in constructor
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)
{
return ( A.size() == 0 );
......@@ -158,6 +152,8 @@ namespace Dune {
return llsMatrixPseudoInv;
}
// TODO: avoid usage of '_' in parameter name - either use
// static method or have no parameters
Matrix matrixSetUp(const Matrix &llsMatrix_, const Matrix &constraintMatrix_)
{
if ( isEmpty(constraintMatrix_) ) {
......@@ -200,8 +196,6 @@ namespace Dune {
// std::cout << "System matrix: " << std::endl;
// printMatrix(systemMatrix);
assert(llsMatrix_.cols() + 1 < systemMatrix.size());
systemMatrix.invert();
// std::cout << "System matrix invert: " << std::endl;
......@@ -223,7 +217,7 @@ namespace Dune {
llsMatrix_.usmtv(2,b,y);
for (Size i = 0; i < systemVector.size(); ++i){
if (i < llsMatrix_.cols() ) {
if (i < y.size() ) {
systemVector[i] = y[i];
}
else {
......
......@@ -287,13 +287,14 @@ namespace Dune
assert( v.size() == size() );
for ( std::size_t i=0; i < size(); ++i)
matrix_[i][col_] = v[i];
return *this;
}
Matrix &matrix_;
int col_;
};
template <class Matrix>
ColumnVector<Matrix> columnVector(Matrix &matrix, int col)
{ return ColumnVector(matrix,col); }
{ return ColumnVector<Matrix>(matrix,col); }
// L.solve(d,b,columnVector(valueProjection,beta));
......@@ -384,33 +385,30 @@ namespace Dune
<< 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;
// 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 );
LeftPseudoInverse< DomainFieldType > pseudoInverse( numShapeFunctions );
// these are the matrices we need to compute
valueProjections_.resize( agglomeration().size() );
jacobianProjections_.resize( agglomeration().size() );
hessianProjections_.resize( agglomeration().size() );
stabilizations_.resize( agglomeration().size() );
// 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 )
// 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 );
std::cout << "Num dofs: "
<< numDofs << std::endl;
const auto &bbox = blockMapper_.indexSet().boundingBox(agglomerate);
const std::size_t numDofs = blockMapper().numDofs( agglomerate );
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) );
......@@ -485,43 +483,22 @@ namespace Dune
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 colValueProjection = ColumnVector( valueProjection, beta );
for ( std::size_t beta = 0; beta < numDofs; ++beta )
{
auto colValueProjection = columnVector( valueProjection, beta );
// set up vectors b and d needed for least squares
b[ beta ] = 1;
if( beta >= numDofs - numInnerShapeFunctions ){
d[ beta - numDofs + numInnerShapeFunctions] = 1;
}
std::cout << "b: " << std::endl;
printVector(b);
std::cout << "d: " << std::endl;
printVector(d);
if( beta >= numDofs - numInnerShapeFunctions )
d[ beta - numDofs + numInnerShapeFunctions] = H0;
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;
if( beta >= numDofs - numInnerShapeFunctions )
d[ beta - numDofs + numInnerShapeFunctions] = 0;
b[ beta ] = 0;
}
#if 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) {
......@@ -529,8 +506,9 @@ namespace Dune
}
std::cout << std::endl;
}
#endif
std::cout << "I reached here 2" << std::endl;
/////////////////////////////////////////
/////////////////////////////////////////
// !!! Original value projection implementation
......
......@@ -26,14 +26,14 @@ useGrid = 0 # 0..3
# Note: suboptimal laplace error for bubble (space is reduced to polorder=3 but could be 4 = ts+2
methods = [ ### "[space,scheme,spaceKwrags]"
["lagrange","galerkin",{}, "Lagrange"],
["vem","vem",{"testSpaces":[ [0], [order-2], [order-1] ] }, "Bubble"],
["vem","vem",{"testSpaces":[ [0], [order-2], [order-3] ] }, "Serendipity"],
# ["vem","vem",{"testSpaces":[ [0], [order-2], [order-1] ] }, "Bubble"],
# ["vem","vem",{"testSpaces":[ [0], [order-2], [order-3] ] }, "Serendipity"],
# ["vem","vem",{"testSpaces":[ [-1], [order-1], [order-3] ] }, "Nc-Serendipity"],
["vem","vem",{"conforming":True}, "conforming"],
["vem","vem",{"conforming":False}, "non-conforming"],
["vem","vem",{"testSpaces":[ [0], [order-3,order-2], [order-4] ] }, "C1-non-conforming"],
["vem","vem",{"testSpaces":[ [0], [order-2,order-2], [order-2] ] }, "C1C0-conforming"],
["bbdg","bbdg",{},"bbdg"],
# ["vem","vem",{"testSpaces":[ [0], [order-3,order-2], [order-4] ] }, "C1-non-conforming"],
# ["vem","vem",{"testSpaces":[ [0], [order-2,order-2], [order-2] ] }, "C1C0-conforming"],
# ["bbdg","bbdg",{},"bbdg"],
]
uflSpace = dune.ufl.Space(2, dimRange=1)
......@@ -91,7 +91,9 @@ for level in range(2,5):
figPos = 100*figRows+10*figCols+1
res = []
for i,m in enumerate(methods):
dune.generator.setFlags(flags="-g",noChecks=True)
space = create.space(m[0], polyGrid, order=order, dimRange=1, storage="istl", **m[2])
dune.generator.unsetFlags()
dfs,errors,info = compute(polyGrid, space, m[1])
print("method:(",m[0],m[2],")",
"Size: ",space.size,
......
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