Skip to content
Snippets Groups Projects
Commit 59cb3c3d authored by Carsten Gräser's avatar Carsten Gräser
Browse files

[examples] Update Hermite examples

parent 8d4a0883
No related branches found
No related tags found
No related merge requests found
Pipeline #77193 canceled
Stage: .pre
Stage: test
Hide whitespace changes
Inline Side-by-side
examples/biharmonic.cc
+ 29
− 18
View file @ 59cb3c3d
......@@ -40,6 +40,9 @@
#include "utilities.hh"
// A custom differential operator for evaluating the Hessian
// that simply uses an internal vector to cache the Hessian
// values at quadrature points.
template<class B, class TP, std::size_t argIndex>
class FEFunctionHessianOperator : public Dune::Fufem::Forms::FEOperatorBase<B, TP, argIndex>
{
......@@ -149,7 +152,6 @@ int main (int argc, char *argv[]) try
grid->globalRefine(refinements);
auto gridView = grid->leafGridView();
using GridView = decltype(gridView);
/////////////////////////////////////////////////////////
// Choose a finite element space
......@@ -164,7 +166,7 @@ int main (int argc, char *argv[]) try
/////////////////////////////////////////////////////////
using Vector = Dune::BlockVector<double>;
using BitVector = Dune::BlockVector<double>;
using BitVector = std::vector<bool>;
using Matrix = Dune::BCRSMatrix<double>;
Vector rhs;
......@@ -193,23 +195,30 @@ int main (int argc, char *argv[]) try
// return 0.0;
// };
auto exactSol = [&](const auto& x) { return 16*x[0]*(x[0]-1)*x[1]*(x[1]-1); };
auto rightHandSide = [] (const auto& x) { return 8*16;};
auto dirichletValuesC0 = [&](const auto& x){
return exactSol(x);
};
// auto exactSol = [&](const auto& x) { return 16*x[0]*(x[0]-1)*x[1]*(x[1]-1); };
// auto rightHandSide = [] (const auto& x) { return 8*16;};
// auto dirichletValuesC0 = [&](const auto& x){
// return exactSol(x);
// };
// auto exactSol = [&](const auto& x) { return 0; };
// auto rightHandSide = [] (const auto& x) { return 10;};
// auto dirichletValuesC0 = [pi](const auto& x){
// using std::sin;
// return sin(2*pi*x[0]);
// };
auto rightHandSide = [] (const auto& x) {
return 500;
};
auto dirichletValuesC0 = [pi](const auto& x){
return 0.0;
};
// Make Dirichlet data differentiable, since this is required by local interpolation
using SignatureTag = Dune::Functions::SignatureTag<double(Dune::FieldVector<double,2>)>;
auto dirichletValues = Dune::Fufem::makeAdolCFunction(SignatureTag(), dirichletValuesC0);
// Disable parallel assembly for UGGrid before 2.10
// Disable parallel assembly for UGGrid before 2.10
#if DUNE_VERSION_GT(DUNE_FUNCTIONS, 2, 9)
auto gridViewPartition = Dune::Fufem::coloredGridViewPartition(gridView);
auto globalAssembler = GlobalAssembler(basis, basis, 1, std::cref(gridViewPartition), threadCount);
......@@ -225,7 +234,6 @@ int main (int argc, char *argv[]) try
auto u = trialFunction(basis, NonAffineFamiliy{});
auto f = Coefficient(Functions::makeGridViewFunction(rightHandSide, gridView));
// auto a = integrate(dot(grad(u), grad(v)));
auto a = integrate(dot(hessian(u), hessian(v)));
auto b = integrate(f*v);
......@@ -247,6 +255,7 @@ int main (int argc, char *argv[]) try
Dune::Functions::istlVectorBackend(isConstrained).resize(basis);
Dune::Functions::istlVectorBackend(isConstrained) = 0;
// Mark all boundary DOFs with derivative order one in normal direction
{
auto normalDerivativeOrder = [](const auto& functionalDescriptor, const auto& unitNormal) -> std::size_t
{
......@@ -283,8 +292,7 @@ int main (int argc, char *argv[]) try
interpolate(basis, sol, dirichletValues, isConstrained);
incorporateEssentialConstraints(basis, stiffnessMatrix, rhs, isConstrained, sol);
incorporateEssentialConstraints(stiffnessMatrix, rhs, isConstrained, sol);
////////////////////////////
// Compute solution
......@@ -293,14 +301,19 @@ int main (int argc, char *argv[]) try
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> op(stiffnessMatrix);
// FastAMG is not working for non-blocked matrices in 2.10
#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 11)
auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
#else
// Sequential incomplete LU decomposition as the preconditioner
SeqILDL<Matrix,Vector,Vector> ildl(stiffnessMatrix,1.0);
auto preconditioner = SeqILDL<Matrix,Vector,Vector>(stiffnessMatrix,1.0);
#endif
// Preconditioned conjugate-gradient solver
CGSolver<Vector> cg(op,
ildl, // preconditioner
preconditioner, // preconditioner
1e-4, // desired residual reduction factor
100, // maximum number of iterations
1000, // maximum number of iterations
2); // verbosity of the solver
// Object storing some statistics about the solving process
......@@ -314,7 +327,6 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////////////////////////////////////
auto solFunction = Functions::makeDiscreteGlobalBasisFunction<double>(basis, sol);
auto exactSolFunction = Functions::makeGridViewFunction(exactSol, gridView);
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
......@@ -329,15 +341,14 @@ int main (int argc, char *argv[]) try
auto vtkWriter = UnstructuredGridWriter(DiscontinuousLagrangeDataCollector(basis.gridView(), 3));
vtkWriter.addPointData(u, VTK::FieldInfo("sol", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.addPointData(jacobian(u), VTK::FieldInfo("grad_sol", VTK::FieldInfo::Type::vector, 2));
vtkWriter.addPointData(exactSolFunction, VTK::FieldInfo("sol_exact", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("biharmonic");
#else
// We need to use the SubsamplingVTKWriter, because dune-vtk's LagrangeDataCollector
// does not work with mixed grids in 2.9.
using GridView = decltype(gridView);
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));
vtkWriter.addVertexData(u, VTK::FieldInfo("sol", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.addVertexData(grad(u), VTK::FieldInfo("grad_sol", VTK::FieldInfo::Type::vector, 2));
vtkWriter.addVertexData(exactSolFunction, VTK::FieldInfo("sol_exact", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("biharmonic");
#endif
......
examples/poisson-hermite-nitsche.cc
+ 2
− 7
View file @ 59cb3c3d
......@@ -20,9 +20,8 @@
#include <dune/istl/solvers.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/cubichermitebasis.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/facenormalgridfunction.hh>
......@@ -35,8 +34,6 @@
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#endif
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/fufem/parallel/elementcoloring.hh>
#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
#include <dune/fufem/assemblers/dunefunctionsfunctionalassembler.hh>
......@@ -216,8 +213,6 @@ int main (int argc, char *argv[]) try
std::cout << "Number of DOFs is " << feBasis.dimension() << std::endl;
auto dirichletIndicatorFunction = [](auto x) { return true; };
auto boundary = BoundaryPatch(feBasis.gridView(), true);
auto rightHandSide = [] (const auto& x) { return 10;};
......@@ -225,7 +220,7 @@ int main (int argc, char *argv[]) try
auto dirichletValues = [pi](const auto& x){ return std::sin(2*pi*x[0]); };
// Disable parallel assembly for UGGrid before 2.10
// Disable parallel assembly for UGGrid before 2.10
#if DUNE_VERSION_GT(DUNE_FUNCTIONS, 2, 9)
auto gridViewPartition = Dune::Fufem::coloredGridViewPartition(gridView);
auto globalAssembler = GlobalAssembler(feBasis, feBasis, 1, std::cref(gridViewPartition), threadCount);
......
examples/poisson-hermite.cc
+ 42
− 13
View file @ 59cb3c3d
......@@ -32,6 +32,7 @@
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#endif
#include <dune/fufem/logger.hh>
#include <dune/fufem/functions/adolcfunction.hh>
#include <dune/fufem/parallel/elementcoloring.hh>
#include <dune/fufem/boundarypatch.hh>
......@@ -46,6 +47,9 @@ using namespace Dune;
int main (int argc, char *argv[]) try
{
auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
std::size_t threadCount = 4;
if (argc>1)
threadCount = std::stoul(std::string(argv[1]));
......@@ -54,18 +58,24 @@ int main (int argc, char *argv[]) try
if (argc>2)
refinements = std::stoul(std::string(argv[2]));
log("Command line processed");
// Set up MPI, if available
MPIHelper::instance(argc, argv);
log("MPI initialized");
///////////////////////////////////
// Generate the grid
///////////////////////////////////
auto grid = Dune::StructuredGridFactory<Dune::UGGrid<2>>::createSimplexGrid({0,0}, {1, 1}, {2, 2});
log("Grid created");
grid->globalRefine(refinements);
log(Dune::formatString("Grid refined %d times", refinements));
auto gridView = grid->leafGridView();
using GridView = decltype(gridView);
......@@ -75,15 +85,17 @@ int main (int argc, char *argv[]) try
using namespace Functions::BasisFactory;
// auto basis = makeBasis(gridView, cubicHermite());
auto basis = makeBasis(gridView, reducedCubicHermite());
auto basis = makeBasis(gridView, cubicHermite());
// auto basis = makeBasis(gridView, reducedCubicHermite());
log(Dune::formatString("Basis created with dimension %d", basis.dimension()));
/////////////////////////////////////////////////////////
// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
using Vector = Dune::BlockVector<double>;
using BitVector = Dune::BlockVector<double>;
using BitVector = std::vector<bool>;
using Matrix = Dune::BCRSMatrix<double>;
Vector rhs;
......@@ -95,8 +107,6 @@ int main (int argc, char *argv[]) try
// Assemble the system
/////////////////////////////////////////////////////////
std::cout << "Number of DOFs is " << basis.dimension() << std::endl;
auto dirichletIndicatorFunction = [](auto x) { return true; };
auto dirichletPatch = BoundaryPatch(basis.gridView());
......@@ -109,15 +119,18 @@ int main (int argc, char *argv[]) try
return sin(2*pi*x[0]);
};
// Make Dirichlet data differentiable, since this is required by local interpolation
using SignatureTag = Dune::Functions::SignatureTag<double(Dune::FieldVector<double,2>)>;
auto dirichletValues = Dune::Fufem::makeAdolCFunction(SignatureTag(), dirichletValuesC0);
// Disable parallel assembly for UGGrid before 2.10
// Disable parallel assembly for UGGrid before 2.10
#if DUNE_VERSION_GT(DUNE_FUNCTIONS, 2, 9)
auto gridViewPartition = Dune::Fufem::coloredGridViewPartition(gridView);
auto globalAssembler = GlobalAssembler(basis, basis, 1, std::cref(gridViewPartition), threadCount);
log("Grid coloring computed");
auto globalAssembler = GlobalAssembler(basis, basis, 0, std::cref(gridViewPartition), threadCount);
#else
auto globalAssembler = GlobalAssembler(basis, basis, 1);
auto globalAssembler = GlobalAssembler(basis, basis, 0);
#endif
{
......@@ -131,8 +144,11 @@ int main (int argc, char *argv[]) try
auto a = integrate(dot(grad(u), grad(v)));
auto b = integrate(f*v);
log("Assembler set up");
globalAssembler.assembleOperator(stiffnessMatrix, a);
log("Matrix assembled");
globalAssembler.assembleFunctional(rhs, b);
log("RHS assembled");
}
// *********************************************
......@@ -149,6 +165,7 @@ int main (int argc, char *argv[]) try
Dune::Functions::istlVectorBackend(isConstrained).resize(basis);
Dune::Functions::istlVectorBackend(isConstrained) = 0;
// Mark all boundary DOFs with derivative order zero in normal direction
{
auto normalDerivativeOrder = [](const auto& functionalDescriptor, const auto& unitNormal) -> std::size_t
{
......@@ -184,9 +201,10 @@ int main (int argc, char *argv[]) try
}
interpolate(basis, sol, dirichletValues, isConstrained);
log("Boundary DOFs interpolated");
incorporateEssentialConstraints(basis, stiffnessMatrix, rhs, isConstrained, sol);
incorporateEssentialConstraints(stiffnessMatrix, rhs, isConstrained, sol);
log("Boundary condition incorporated into system");
////////////////////////////
// Compute solution
......@@ -195,22 +213,33 @@ int main (int argc, char *argv[]) try
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> op(stiffnessMatrix);
// FastAMG is not working for non-blocked matrices in 2.10
#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 11)
auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
#else
// Sequential incomplete LU decomposition as the preconditioner
SeqILDL<Matrix,Vector,Vector> ildl(stiffnessMatrix,1.0);
auto preconditioner = SeqILDL<Matrix,Vector,Vector>(stiffnessMatrix,1.0);
#endif
log("Preconditioner created");
// Preconditioned conjugate-gradient solver
CGSolver<Vector> cg(op,
ildl, // preconditioner
preconditioner, // preconditioner
1e-4, // desired residual reduction factor
100, // maximum number of iterations
2); // verbosity of the solver
log("Solver created");
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
// Solve!
cg.apply(sol, rhs, statistics);
log("Linear system solved");
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
////////////////////////////////////////////////////////////////////////////
......
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