From 59cb3c3da25ebb8129f66d0d0b46f69893f7bf0d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@math.fau.de>
Date: Mon, 3 Mar 2025 14:55:16 +0100
Subject: [PATCH] [examples] Update Hermite examples
---
examples/biharmonic.cc | 47 ++++++++++++++----------
examples/poisson-hermite-nitsche.cc | 9 ++---
examples/poisson-hermite.cc | 55 ++++++++++++++++++++++-------
3 files changed, 73 insertions(+), 38 deletions(-)
diff --git a/examples/biharmonic.cc b/examples/biharmonic.cc
index 89584d74..bb122c7e 100644
--- a/examples/biharmonic.cc
+++ b/examples/biharmonic.cc
@@ -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
diff --git a/examples/poisson-hermite-nitsche.cc b/examples/poisson-hermite-nitsche.cc
index cd18d5a4..0cc871a6 100644
--- a/examples/poisson-hermite-nitsche.cc
+++ b/examples/poisson-hermite-nitsche.cc
@@ -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);
diff --git a/examples/poisson-hermite.cc b/examples/poisson-hermite.cc
index 14ddc46b..af9ef637 100644
--- a/examples/poisson-hermite.cc
+++ b/examples/poisson-hermite.cc
@@ -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
////////////////////////////////////////////////////////////////////////////
--
GitLab