Skip to content
Snippets Groups Projects
Commit 64d05ceb authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Assemble into a MultiTypeBlockMatrix 

The stiffness matrix has a 2x2 block structure, and three of the
four blocks are dense.  Use a MultiTypeBlockMatrix data type with
sparse and dense subblocks to properly represent this structure.
parent 09e83287
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/microstructure/CorrectorComputer.hh
+ 78
− 74
View file @ 64d05ceb
......@@ -5,7 +5,12 @@
#include <dune/common/parametertree.hh>
// #include <dune/common/float_cmp.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
......@@ -38,8 +43,19 @@ public:
using FuncVector = std::function< VectorRT(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using Func2int = std::function< int(const Domain&) >;
using VectorCT = Dune::BlockVector<Dune::FieldVector<double,1> >;
using MatrixCT = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
using V0 = Dune::BlockVector<double>;
using V1 = Dune::BlockVector<double>;
using VectorCT = Dune::MultiTypeBlockVector<V0, V1>;
using SparseMatrix = Dune::BCRSMatrix<double>;
using DenseMatrix = Dune::Matrix<double>;
using MatrixCT = Dune::MultiTypeBlockMatrix<
Dune::MultiTypeBlockVector<SparseMatrix, DenseMatrix>,
Dune::MultiTypeBlockVector<DenseMatrix, DenseMatrix>
>;
using ElementMatrixCT = Dune::Matrix<double>;
protected:
......@@ -61,7 +77,7 @@ protected:
VectorCT load_alpha1_,load_alpha2_,load_alpha3_; //right-hand side(load) vectors
VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
VectorCT phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
Dune::BlockVector<double> phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
Dune::FieldVector<double,3> m_1_, m_2_, m_3_; // Corrector m_i coefficient vectors
// (assembled) corrector matrices M_i
......@@ -87,9 +103,6 @@ protected:
const std::array<Func2Tensor, 3> x3MatrixBasisContainer_ = {x3G_1_, x3G_2_, x3G_3_};
// --- Offset between basis indices
const int phiOffset_;
public:
///////////////////////////////
// constructor
......@@ -106,8 +119,7 @@ public:
parameterSet_(parameterSet),
matrixBasis_(std::array<VoigtVector<double,3>,3>{matrixToVoigt(Dune::FieldMatrix<double,3,3>({{1, 0, 0}, {0, 0, 0}, {0, 0, 0}})),
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))}),
phiOffset_(basis.size())
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))})
{
assemble();
......@@ -166,9 +178,9 @@ public:
// shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
auto getCorr_phi1(){return std::make_shared<Dune::BlockVector<double>>(phi_1_);}
auto getCorr_phi2(){return std::make_shared<Dune::BlockVector<double>>(phi_2_);}
auto getCorr_phi3(){return std::make_shared<Dune::BlockVector<double>>(phi_3_);}
......@@ -180,9 +192,8 @@ public:
// | phi*phi m*phi |
// | phi *m m*m |
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
nb.resize(basis_.size()+3, basis_.size()+3);
nb.resize(basis_.size(), basis_.size());
for (const auto& element : elements(basis_.gridView()))
{
......@@ -198,21 +209,6 @@ public:
nb.add(row[0],col[0]); // nun würde auch nb.add(row,col) gehen..
}
}
for (size_t i=0; i<localView.size(); i++)
{
auto row = localView.index(i);
for (size_t j=0; j<3; j++ )
{
nb.add(row,phiOffset+j); // m*phi part of StiffnessMatrix
nb.add(phiOffset+j,row); // phi*m part of StiffnessMatrix
}
}
for (size_t i=0; i<3; i++ )
for (size_t j=0; j<3; j++ )
{
nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
}
}
//////////////////////////////////////////////////////////////////
// setOneBaseFunctionToZero
......@@ -627,13 +623,21 @@ public:
{
std::cout << "assemble Stiffness-Matrix begins." << std::endl;
using namespace Dune::Indices; // For _0, _1, etc.
// Set up the sparse block
Dune::MatrixIndexSet occupationPattern;
getOccupationPattern(occupationPattern);
occupationPattern.exportIdx(matrix);
occupationPattern.exportIdx(matrix[_0][_0]);
// Set up the three dense blocks
matrix[_0][_1].setSize(basis_.size(), 3);
matrix[_1][_0].setSize(3, basis_.size());
matrix[_1][_1].setSize(3,3);
matrix = 0;
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
// A cache for the element stiffness matrices, if desired
bool cacheElementMatrices = true;
......@@ -721,18 +725,18 @@ public:
{
auto row = localView.index(i);
auto col = localView.index(j);
matrix[row][col] += elementMatrix[i][j];
matrix[_0][_0][row][col] += elementMatrix[i][j];
}
for (size_t i=0; i<localPhiOffset; i++)
for(size_t m=0; m<3; m++)
{
auto row = localView.index(i);
matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m];
matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
matrix[_0][_1][row][m] += elementMatrix[i][localPhiOffset+m];
matrix[_1][_0][m][row] += elementMatrix[localPhiOffset+m][i];
}
for (size_t m=0; m<3; m++ )
for (size_t n=0; n<3; n++ )
matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
matrix[_1][_1][m][n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
// printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
}
......@@ -744,11 +748,12 @@ public:
)
{
// std::cout << "assemble load vector." << std::endl;
b.resize(basis_.size()+3);
using namespace Dune::Indices;
b[_0].resize(basis_.size());
b[_1].resize(3);
b = 0;
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
// Transform G_alpha's to GridViewFunctions/LocalFunctions
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis_.gridView());
......@@ -771,7 +776,7 @@ public:
const auto localPhiOffset = localView.size();
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
VectorCT elementRhs;
Dune::BlockVector<double> elementRhs;
// std::cout << "----------------------------------Element : " << counter << std::endl; //TEST
// counter++;
// computeElementLoadVector(localView, muLocal, lambdaLocal, elementRhs, loadFunctional);
......@@ -785,10 +790,10 @@ public:
for (size_t p=0; p<localPhiOffset; p++)
{
auto row = localView.index(p);
b[row] += elementRhs[p];
b[_0][row] += elementRhs[p];
}
for (size_t m=0; m<3; m++ )
b[phiOffset+m] += elementRhs[localPhiOffset+m];
b[_1][m] += elementRhs[localPhiOffset+m];
//printvector(std::cout, b, "b", "--");
}
}
......@@ -837,15 +842,20 @@ public:
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
Dune::FieldVector<std::size_t,3> row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
using namespace Dune::Indices; // For _0, _1, etc.
for (int k = 0; k<dim; k++)
{
load_alpha1_[row[k]] = 0.0;
load_alpha2_[row[k]] = 0.0;
load_alpha3_[row[k]] = 0.0;
auto cIt = stiffnessMatrix_[row[k]].begin();
auto cEndIt = stiffnessMatrix_[row[k]].end();
load_alpha1_[_0][row[k]] = 0.0;
load_alpha2_[_0][row[k]] = 0.0;
load_alpha3_[_0][row[k]] = 0.0;
auto cIt = stiffnessMatrix_[_0][_0][row[k]].begin();
auto cEndIt = stiffnessMatrix_[_0][_0][row[k]].end();
for (; cIt!=cEndIt; ++cIt)
*cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
// Delete corresponding row in off-diagonal block
stiffnessMatrix_[_0][_1][row[k]] = 0.0;
}
}
......@@ -886,7 +896,7 @@ public:
}
auto integralMean(VectorCT& coeffVector)
auto integralMean(const Dune::BlockVector<double>& coeffVector)
{
auto GVFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,dim>>(basis_,coeffVector);
auto localfun = localFunction(GVFunction);
......@@ -921,7 +931,7 @@ public:
}
auto subtractIntegralMean(VectorCT& coeffVector)
auto subtractIntegralMean(Dune::BlockVector<double>& coeffVector)
{
// Subtract correct integral mean from each associated component function
auto IM = integralMean(coeffVector);
......@@ -932,6 +942,7 @@ public:
auto idx = childToIndexMap(k);
for ( int i : idx)
coeffVector[i] -= IM[k];
// TODO: Do we have to modify coeffVector[_1] as well?
}
}
......@@ -1100,25 +1111,16 @@ public:
////////////////////////////////////////////////////////////////////////////////////
// Extract phi_alpha & M_alpha coefficients
////////////////////////////////////////////////////////////////////////////////////
phi_1_.resize(basis_.size());
phi_1_ = 0;
phi_2_.resize(basis_.size());
phi_2_ = 0;
phi_3_.resize(basis_.size());
phi_3_ = 0;
for(size_t i=0; i<basis_.size(); i++)
{
phi_1_[i] = x_1_[i];
phi_2_[i] = x_2_[i];
phi_3_[i] = x_3_[i];
}
for(size_t i=0; i<3; i++)
{
m_1_[i] = x_1_[phiOffset_+i];
m_2_[i] = x_2_[phiOffset_+i];
m_3_[i] = x_3_[phiOffset_+i];
}
using namespace Dune::Indices;
phi_1_ = x_1_[_0];
phi_2_ = x_2_[_0];
phi_3_ = x_3_[_0];
std::copy(x_1_[_1].begin(), x_1_[_1].end(), m_1_.begin());
std::copy(x_2_[_1].begin(), x_2_[_1].end(), m_2_.begin());
std::copy(x_3_[_1].begin(), x_3_[_1].end(), m_3_.begin());
// assemble M_alpha's (actually iota(M_alpha) )
for (auto& matrix : mContainer)
......@@ -1160,9 +1162,9 @@ public:
subtractIntegralMean(phi_1_);
subtractIntegralMean(phi_2_);
subtractIntegralMean(phi_3_);
subtractIntegralMean(x_1_);
subtractIntegralMean(x_2_);
subtractIntegralMean(x_3_);
subtractIntegralMean(x_1_[_0]);
subtractIntegralMean(x_2_[_0]);
subtractIntegralMean(x_3_[_0]);
//////////////////////////////////////////
// Check Integral-mean again:
//////////////////////////////////////////
......@@ -1182,19 +1184,19 @@ public:
{
log_ << "\nSolution of Corrector problems:\n";
log_ << "\n Corrector_phi1:\n";
log_ << x_1_ << std::endl;
log_ << x_1_[_0] << std::endl;
}
if(parameterSet_.get<bool>("write_corrector_phi2", false))
{
log_ << "-----------------------------------------------------";
log_ << "\n Corrector_phi2:\n";
log_ << x_2_ << std::endl;
log_ << x_2_[_0] << std::endl;
}
if(parameterSet_.get<bool>("write_corrector_phi3", false))
{
log_ << "-----------------------------------------------------";
log_ << "\n Corrector_phi3:\n";
log_ << x_3_ << std::endl;
log_ << x_3_[_0] << std::endl;
}
}
......@@ -1443,6 +1445,8 @@ public:
auto localView = basis_.localView();
using namespace Dune::Indices; // For _0, _1, etc.
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
......@@ -1453,21 +1457,21 @@ public:
{
auto row = localView.index(i);
auto col = localView.index(j);
if(abs( stiffnessMatrix_[row][col] - stiffnessMatrix_[col][row]) > 1e-12 )
if(abs( stiffnessMatrix_[_0][_0][row][col] - stiffnessMatrix_[_0][_0][col][row]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
for(size_t i=0; i<localPhiOffset; i++)
for(size_t m=0; m<3; m++)
{
auto row = localView.index(i);
if(abs( stiffnessMatrix_[row][phiOffset_+m] - stiffnessMatrix_[phiOffset_+m][row]) > 1e-12 )
if(abs( stiffnessMatrix_[_0][_1][row][m] - stiffnessMatrix_[_1][_0][m][row]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
for(size_t m=0; m<3; m++ )
for(size_t n=0; n<3; n++ )
{
if(abs(stiffnessMatrix_[phiOffset_+m][phiOffset_+n] - stiffnessMatrix_[phiOffset_+n][phiOffset_+m]) > 1e-12 )
if(abs(stiffnessMatrix_[_1][_1][m][n] - stiffnessMatrix_[_1][_1][n][m]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
}
......
dune/microstructure/EffectiveQuantitiesComputer.hh
+ 1
− 1
View file @ 64d05ceb
......@@ -91,7 +91,7 @@ public:
auto basis = *correctorComputer_.getBasis();
Dune::ParameterTree parameterSet = correctorComputer_.getParameterSet();
shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(),
shared_ptr<Dune::BlockVector<double>> phiBasis[3] = {correctorComputer_.getCorr_phi1(),
correctorComputer_.getCorr_phi2(),
correctorComputer_.getCorr_phi3()
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment