Commit d4b4cbba authored by Jakub Both's avatar Jakub Both

Change Poisson test case for MFEM.

Use new flux BC to create a nonlinear pressure field.
Allow for higher order RT and BDM elements.

Working:
2D: RT0 x P0, RT1 x P1, BDM1 x P0, BDM2 x P1
3D: RT0 x P0

Not working:
2D: RT2 x P2 - solver does not converge.
3D: RT1 x P1 - fluxes in z direction are not continuous. Rest looks good.
3D: BDM1 x P0 - no interpolation provided by dune-localfunctions.
parent 3156c886
......@@ -28,9 +28,9 @@
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#define RT0
//#define BDM1
//#define DIM2
#define RT
//#define BDM
#define DIM2
using namespace Dune;
......@@ -43,7 +43,7 @@ void getLocalMatrix(const LocalView& localView,
typedef typename LocalView::Element Element;
const Element& element = localView.element();
const int dim = Element::dimension;
const int dim = Element::dimension;
auto geometry = element.geometry();
// Set all matrix entries to zero
......@@ -60,44 +60,66 @@ void getLocalMatrix(const LocalView& localView,
const auto& quad = QuadratureRules<double, dim>::rule(element.type(), order);
// Loop over all quadrature points
for (const auto& quadPoint : quad) {
for (const auto& quadPoint : quad)
{
// Position of the current quadrature point in the reference element
const auto quadPos = quadPoint.position();
// Values of the pressure shape functions
std::vector<FieldVector<double,1> > pressureValues;
pressureLocalFiniteElement.localBasis().evaluateFunction(quadPos, pressureValues);
// The transposed (inverse) Jacobian of the map from the reference element to the element
const auto& jacobianTransposed = geometry.jacobianTransposed(quadPos);
// Values of the flux shape functions
std::vector<FieldVector<double,dim> > fluxValues;
fluxLocalFiniteElement.localBasis().evaluateFunction(quadPos, fluxValues);
// The multiplicative factor in the integral transformation formula and its inverse
const auto integrationElement = geometry.integrationElement(quadPos);
const auto invIntegrationElement = 1.0 / integrationElement;
// The transposed inverse Jacobian of the map from the reference element to the element
const auto& jacobian = geometry.jacobianInverseTransposed(quadPos);
///////////////////////////////////////////////////////////////////////////
// Shape functions - flux - [assume H(div) elements]
///////////////////////////////////////////////////////////////////////////
// The multiplicative factor in the integral transformation formula
const double integrationElement = geometry.integrationElement(quadPos);
// Values of the flux shape functions on the reference element
std::vector<FieldVector<double,dim> > fluxReferenceValues(fluxLocalFiniteElement.size());
fluxLocalFiniteElement.localBasis().evaluateFunction(quadPos, fluxReferenceValues);
// Values of the shape functions on the current element
std::vector<FieldVector<double,dim> > fluxValues(fluxLocalFiniteElement.size());
for (size_t i=0; i<fluxValues.size(); i++)
{
jacobianTransposed.mtv(fluxReferenceValues[i], fluxValues[i]);
fluxValues[i] *= invIntegrationElement;
}
// Gradients of the flux shape functions on the reference element
std::vector<FieldMatrix<double,dim,dim> > fluxReferenceGradients;
std::vector<FieldMatrix<double,dim,dim> > fluxReferenceGradients(fluxLocalFiniteElement.size());
fluxLocalFiniteElement.localBasis().evaluateJacobian(quadPos, fluxReferenceGradients);
// Compute the flux shape function gradients on the real element
std::vector<FieldMatrix<double,dim,dim> > fluxGradients(fluxReferenceGradients.size());
for (size_t i=0; i<fluxGradients.size(); i++)
// Compute the flux shape function divergence on the real element
std::vector<double> fluxDivergence(fluxValues.size(), 0.0);
for (size_t i=0; i<fluxValues.size(); i++)
{
for (size_t j=0; j<dim; j++)
jacobian.mv(fluxReferenceGradients[i][j], fluxGradients[i][j]);
fluxDivergence[i] += fluxReferenceGradients[i][j][j];
fluxDivergence[i] *= invIntegrationElement;
}
///////////////////////////////////////////////////////////////////////////
// Shape functions - pressure
///////////////////////////////////////////////////////////////////////////
// Values of the pressure shape functions
std::vector<FieldVector<double,1> > pressureValues(pressureLocalFiniteElement.size());
pressureLocalFiniteElement.localBasis().evaluateFunction(quadPos, pressureValues);
/////////////////////////////
// Flux--flux coupling
/////////////////////////////
for (size_t i=0; i<fluxLocalFiniteElement.size(); i++) {
for (size_t i=0; i<fluxLocalFiniteElement.size(); i++)
{
size_t row = localView.tree().child(_0).localIndex(i);
for (size_t j=0; j<fluxLocalFiniteElement.size(); j++) {
size_t col = localView.tree().child(_0).localIndex(j);
elementMatrix[row][col] += (fluxValues[i] * fluxValues[j]) * quadPoint.weight() * integrationElement;
for (size_t j=0; j<fluxLocalFiniteElement.size(); j++)
{
size_t col = localView.tree().child(_0).localIndex(j);
elementMatrix[row][col] += (fluxValues[i] * fluxValues[j]) * quadPoint.weight() * integrationElement;
}
}
......@@ -105,19 +127,15 @@ void getLocalMatrix(const LocalView& localView,
// Flux--pressure coupling
/////////////////////////////
for (size_t i=0; i<fluxLocalFiniteElement.size(); i++) {
for (size_t i=0; i<fluxLocalFiniteElement.size(); i++)
{
size_t fluxIndex = localView.tree().child(_0).localIndex(i);
// Pre-compute divergence
double fluxDivergence = 0.;
for (size_t k=0; k<dim; k++)
fluxDivergence += fluxGradients[i][k][k];
for (size_t j=0; j<pressureLocalFiniteElement.size(); j++) {
for (size_t j=0; j<pressureLocalFiniteElement.size(); j++)
{
size_t pressureIndex = localView.tree().child(_1).localIndex(j);
// Pre-compute matrix contribution
double tmp = - (pressureValues[j] * fluxDivergence) * quadPoint.weight() * integrationElement;
double tmp = - (fluxDivergence[i] * pressureValues[j]) * quadPoint.weight() * integrationElement;
elementMatrix[fluxIndex][pressureIndex] += tmp;
elementMatrix[pressureIndex][fluxIndex] += tmp;
......@@ -155,8 +173,8 @@ void getVolumeTerm( const LocalView& localView,
const auto& quad = QuadratureRules<double, dim>::rule(element.type(), order);
// Loop over all quadrature points
for (const auto& quadPoint : quad) {
for (const auto& quadPoint : quad)
{
// Position of the current quadrature point in the reference element
const auto& quadPos = quadPoint.position();
......@@ -170,10 +188,10 @@ void getVolumeTerm( const LocalView& localView,
pressureLocalFiniteElement.localBasis().evaluateFunction(quadPos, pressureValues);
// Actually compute the vector entries
for (size_t j=0; j<pressureLocalFiniteElement.size(); j++) {
for (size_t j=0; j<pressureLocalFiniteElement.size(); j++)
{
size_t pressureIndex = localView.tree().child(_1).localIndex(j);
localRhs[pressureIndex] += (-1) * pressureValues[j] * functionValue * quadPoint.weight() * integrationElement;
localRhs[pressureIndex] += - pressureValues[j] * functionValue * quadPoint.weight() * integrationElement;
}
}
......@@ -202,18 +220,17 @@ void getOccupationPattern(const Basis& basis,
localIndexSet.bind(localView);
// Add element stiffness matrix onto global stiffness matrix
for (size_t i=0; i<localIndexSet.size(); i++) {
for (size_t i=0; i<localIndexSet.size(); i++)
{
// The global index of the i-th local degree of freedom of the element 'e'
auto row = localIndexSet.index(i);
for (size_t j=0; j<localIndexSet.size(); ++j) {
for (size_t j=0; j<localIndexSet.size(); ++j)
{
// The global index set of the j-th local degree of freedom of the element 'e'
auto col = localIndexSet.index(j);
nb[row[0]][col[0]].add(row[1], col[1]);
}
}
......@@ -261,17 +278,16 @@ void assembleMixedPoissonMatrix(const Basis& basis,
getLocalMatrix(localView, elementMatrix);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t i=0; i<elementMatrix.N(); i++) {
for (size_t i=0; i<elementMatrix.N(); i++)
{
// The global index of the i-th local degree of freedom of the element 'e'
auto row = localIndexSet.index(i);
for (size_t j=0; j<elementMatrix.M(); j++ ) {
for (size_t j=0; j<elementMatrix.M(); j++ )
{
// The global index of the j-th local degree of freedom of the element 'e'
auto col = localIndexSet.index(j);
matrix[row[0]][col[0]][row[1]][col[1]] += elementMatrix[i][j];
}
}
......@@ -315,12 +331,11 @@ void assembleMixedPoissonRhs(const Basis& basis,
localVolumeTerm.bind(element);
getVolumeTerm(localView, localRhs, localVolumeTerm);
for (size_t i=0; i<localRhs.size(); i++) {
for (size_t i=0; i<localRhs.size(); i++)
{
// The global index of the i-th vertex of the element 'e'
auto row = localIndexSet.index(i);
rhs[row[0]][row[1]] += localRhs[i];
}
}
......@@ -329,9 +344,9 @@ void assembleMixedPoissonRhs(const Basis& basis,
int main (int argc, char *argv[]) try
{
#ifndef RT0
#ifndef BDM1
DUNE_THROW(Dune::NotImplemented, "Choose RT0 or BDM1.");
#ifndef RT
#ifndef BDM
DUNE_THROW(Dune::NotImplemented, "Choose RT or BDM.");
#endif
#endif
......@@ -344,10 +359,10 @@ int main (int argc, char *argv[]) try
#ifdef DIM2
const int dim = 2;
std::array<int,dim> elements = {20, 20};
std::array<int,dim> elements = {50, 50};
#else
const int dim = 3;
std::array<int,dim> elements = {20, 20, 20};
std::array<int,dim> elements = {30, 30, 30};
#endif
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> l(1);
......@@ -362,17 +377,17 @@ int main (int argc, char *argv[]) try
using namespace Functions::BasisBuilder;
const int k = 0;
auto basis = makeBasis(
gridView,
composite(
#ifdef RT0
rt<0, GeometryType::BasicType::cube>(),
pq<0>()
#endif
#ifdef BDM1
bdm<1, GeometryType::BasicType::cube>(),
pq<0>()
#ifdef RT
rt<k, GeometryType::BasicType::cube>(),
#elif defined BDM
bdm<k+1, GeometryType::BasicType::cube>(),
#endif
pq<k>()
));
......@@ -395,11 +410,12 @@ int main (int argc, char *argv[]) try
assembleMixedPoissonMatrix(basis, stiffnessMatrix);
auto rightHandSide = [] (const Domain& x) { return 2.;};
auto rightHandSide = [] (const Domain& x) { return 10; };
assembleMixedPoissonRhs(basis, rhs, rightHandSide);
auto topFluxBC = [] (const Domain& x) { return 0.0; };
auto lowerFluxBC = [] (const Domain& x) { return 0.0; };
auto pi = std::acos(-1.0);
auto topFluxBC = [&pi] (const Domain& x) { return -0.05 * (1. - x[0]) * std::sin(2.*pi*x[0]); };
auto lowerFluxBC = [&pi] (const Domain& x) { return 0.05 * (1. - x[0]) * std::sin(2.*pi*x[0]); };
using namespace TypeTree::Indices;
using BitVectorType = BlockVector<BlockVector<FieldVector<char,1> > >;
......@@ -408,59 +424,59 @@ int main (int argc, char *argv[]) try
BitVectorType isTopBoundary;
BitVectorType isLowerBoundary;
#ifdef RT0
// Use the same way as for PQk functions to mark boundaries.
// Mark top boundary.
auto topBoundaryIndicator = [&l] (Domain x) {
bool isBoundary = x[dim-1] > l[dim-1] - 1e-8;
return isBoundary;
};
// Mark top boundary.
auto lowerBoundaryIndicator = [&l] (Domain x) {
bool isBoundary = x[dim-1] < 1e-8;
return isBoundary;
};
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicBitVectorView(isTopBoundary), topBoundaryIndicator);
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicBitVectorView(isLowerBoundary), lowerBoundaryIndicator);
#else
// // Use the same way as for PQk functions to mark boundaries - works only for RT0.
//
// // Mark top boundary.
// auto topBoundaryIndicator = [&l] (Domain x)
// {
// bool isBoundary = x[dim-1] > l[dim-1] - 1e-8;
// return isBoundary;
// };
//
// // Mark lower boundary.
// auto lowerBoundaryIndicator = [&l] (Domain x)
// {
// bool isBoundary = x[dim-1] < 1e-8;
// return isBoundary;
// };
//
// interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicBitVectorView(isTopBoundary), topBoundaryIndicator);
// interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicBitVectorView(isLowerBoundary), lowerBoundaryIndicator);
// Use a messy way of defining a boundary. Using the RT0-way, only the RT0 basis DOF are triggered,
// as they suffice to interpolate a constant function on an edge. In particular BDM1 DOF are not triggered,
// as they suffice to interpolate a constant function on an edge. In particular RT1, BDM1 DOF are not triggered,
// which are linear (but non-constant) on edges.
// Mark top boundary.
auto topBoundaryIndicator = [&l] (Domain x) {
double isBoundary = x[dim-1] > l[dim-1] - 1e-8 ? x[0] : 0.0;
return isBoundary;
auto topBoundaryIndicator = [&l] (Domain x)
{
double isBoundary = x[dim-1] > l[dim-1] - 1e-8 ? x[0] : 0.0;
return isBoundary;
};
// Mark top boundary.
auto lowerBoundaryIndicator = [&l] (Domain x) {
double isBoundary = x[dim-1] < 1e-8 ? x[0] : 0.0;
return isBoundary;
// Mark lower boundary.
auto lowerBoundaryIndicator = [&l] (Domain x)
{
double isBoundary = x[dim-1] < 1e-8 ? x[0] : 0.0;
return isBoundary;
};
VectorType isTopBoundaryTmp, isLowerBoundaryTmp;
// Use double-valued interpolation and transfer to char-valued vectors.
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicVectorView(isTopBoundaryTmp), topBoundaryIndicator);
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicVectorView(isTopBoundaryTmp), topBoundaryIndicator);
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicVectorView(isLowerBoundaryTmp), lowerBoundaryIndicator);
HierarchicBitVectorView(isTopBoundary).resize(basis);
HierarchicBitVectorView(isLowerBoundary).resize(basis);
isTopBoundary = 0;
isLowerBoundary = 0;
for (size_t i=0; i<isTopBoundaryTmp[0].size(); i++){
isTopBoundary[0][i] = isTopBoundaryTmp[0][i]!=0 ? 1: 0;
for (size_t i=0; i<isTopBoundaryTmp[0].size(); i++)
{
isTopBoundary[ 0][i] = isTopBoundaryTmp[ 0][i]!=0 ? 1: 0;
isLowerBoundary[0][i] = isLowerBoundaryTmp[0][i]!=0 ? 1: 0;
}
#endif
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicVectorView(rhs), topFluxBC, HierarchicBitVectorView(isTopBoundary));
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicVectorView(rhs), topFluxBC, HierarchicBitVectorView(isTopBoundary));
interpolate(basis, Dune::TypeTree::hybridTreePath(_0), HierarchicVectorView(rhs), lowerFluxBC, HierarchicBitVectorView(isLowerBoundary));
////////////////////////////////////////////
......@@ -468,22 +484,20 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////
// loop over the matrix rows
for (size_t i=0; i<stiffnessMatrix[0][0].N(); i++) {
if (isTopBoundary[0][i] or isLowerBoundary[0][i]) {
for (size_t i=0; i<stiffnessMatrix[0][0].N(); i++)
{
if (isTopBoundary[0][i] or isLowerBoundary[0][i])
{
// Lower right matrix block
auto cIt = stiffnessMatrix[0][0][i].begin();
auto cEndIt = stiffnessMatrix[0][0][i].end();
// loop over nonzero matrix entries in current row
for (; cIt!=cEndIt; ++cIt){
for (; cIt!=cEndIt; ++cIt)
*cIt = (i==cIt.index()) ? 1. : 0.;
}
// Lower left matrix block
for (auto&& entry: stiffnessMatrix[0][1][i])
entry = 0.0;
}
}
......@@ -501,13 +515,9 @@ int main (int argc, char *argv[]) try
// Fancy (but only) way to not have a preconditioner at all
Richardson<VectorType,VectorType> preconditioner(1.0);
// Preconditioned conjugate-gradient solver
RestartedGMResSolver<VectorType> solver(op,
preconditioner,
1e-6, // desired residual reduction factor
100, // number of iterations between restarts
1000, // maximum number of iterations
2); // verbosity of the solver
// Preconditioned GMRES / BiCGSTAB solver
//RestartedGMResSolver<VectorType> solver (op, preconditioner, 1e-6, 1000, 10000, 2);
BiCGSTABSolver<VectorType> solver(op, preconditioner, 1e-6, 10000, 2);
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
......@@ -523,18 +533,21 @@ int main (int argc, char *argv[]) try
using PressureRange = double;
auto fluxFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FluxRange>(basis, TypeTree::hybridTreePath(_0), HierarchicVectorView(x));
auto pressureFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<PressureRange>(basis, TypeTree::hybridTreePath(_1), HierarchicVectorView(x));
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////
SubsamplingVTKWriter<GridView> vtkWriter(gridView,2);
SubsamplingVTKWriter<GridView> vtkWriter(gridView,0);
vtkWriter.addVertexData(fluxFunction, VTK::FieldInfo("flux", VTK::FieldInfo::Type::vector, dim));
vtkWriter.addCellData(pressureFunction, VTK::FieldInfo("pressure", VTK::FieldInfo::Type::scalar, 1));
if (k==0)
vtkWriter.addCellData(pressureFunction, VTK::FieldInfo("pressure", VTK::FieldInfo::Type::scalar, 1));
else
vtkWriter.addVertexData(pressureFunction, VTK::FieldInfo("pressure", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("poisson-mfem-result");
}
// Error handling
catch (Exception e) {
catch (Exception e)
{
std::cout << e << 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