diff --git a/test/sumfact/poisson/facedir-facemod-variation/backup_poisson_dg_unstructured_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/backup_poisson_dg_unstructured_3d_facedir_facemod_variation_driver.cc deleted file mode 100644 index ec5e1ff4646576dd2305e1408bc18dab99ed87cb..0000000000000000000000000000000000000000 --- a/test/sumfact/poisson/facedir-facemod-variation/backup_poisson_dg_unstructured_3d_facedir_facemod_variation_driver.cc +++ /dev/null @@ -1,196 +0,0 @@ -#include "config.h" -#include "dune/common/parallel/mpihelper.hh" -#include "dune/pdelab/stationary/linearproblem.hh" -#include "dune/pdelab/backend/istl.hh" -#include "dune/grid/uggrid.hh" -#include "dune/grid/yaspgrid.hh" -#include "dune/pdelab/finiteelementmap/qkdg.hh" -#include "poisson_dg_unstructured_3d_facedir_facemod_variation_deg1_symdiff_nonquadvec_nongradvec_rOperator_file.hh" -#include "dune/pdelab/gridoperator/gridoperator.hh" -#include "dune/testtools/gridconstruction.hh" -#include "dune/common/parametertree.hh" -#include "dune/common/parametertreeparser.hh" -#include "dune/consistent-edge-orientation/createconsistentgrid.hh" -#include <random> -#include "dune/pdelab/gridfunctionspace/vtk.hh" -#include "dune/grid/io/file/vtk/subsamplingvtkwriter.hh" -#include "string" -#include "dune/codegen/vtkpredicate.hh" -#include "dune/pdelab/function/callableadapter.hh" -#include "dune/pdelab/gridfunctionspace/gridfunctionadapter.hh" -#include "dune/pdelab/common/functionutilities.hh" - - - -int main(int argc, char** argv){ - try - { - // Initialize basic stuff... - Dune::MPIHelper& mpihelper = Dune::MPIHelper::instance(argc, argv); - using RangeType = double; - Dune::ParameterTree initree; - Dune::ParameterTreeParser::readINITree(argv[1], initree); - - // Setup grid (view)... - using Grid = Dune::UGGrid<3>; - - using GV = Grid::LeafGridView; - using DF = Grid::ctype; - IniGridFactory<Grid> factory(initree); - std::shared_ptr<Grid> grid_nonconsistent = factory.getGrid(); - std::shared_ptr<Grid> grid = createConsistentGrid(grid_nonconsistent); - GV gv = grid->leafGridView(); - - // Set up finite element maps... - using DG1_FEM = Dune::PDELab::QkDGLocalFiniteElementMap<DF, RangeType, 1, 3>; - DG1_FEM dg1_fem; - - // Set up grid function spaces... - using VectorBackendDG1 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none>; - using NoConstraintsAssembler = Dune::PDELab::NoConstraints; - using DG1_GFS = Dune::PDELab::GridFunctionSpace<GV, DG1_FEM, NoConstraintsAssembler, VectorBackendDG1>; - DG1_GFS dg1_gfs_(gv, dg1_fem); - dg1_gfs_.name("dg1_gfs_"); - - // Set up constraints container... - using DG1_GFS_CC = DG1_GFS::ConstraintsContainer<RangeType>::Type; - DG1_GFS_CC dg1_gfs__cc; - dg1_gfs__cc.clear(); - Dune::PDELab::constraints(dg1_gfs_, dg1_gfs__cc); - - // Set up grid grid operators... - using LOP_R = rOperator<DG1_GFS, DG1_GFS, RangeType>; - using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>; - using GO_r = Dune::PDELab::GridOperator<DG1_GFS, DG1_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, DG1_GFS_CC, DG1_GFS_CC>; - LOP_R lop_r(dg1_gfs_, dg1_gfs_, initree); - dg1_gfs_.update(); - int generic_dof_estimate = 8 * dg1_gfs_.maxLocalSize(); - int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate); - MatrixBackend mb(dofestimate); - GO_r go_r(dg1_gfs_, dg1_gfs__cc, dg1_gfs_, dg1_gfs__cc, lop_r, mb); - std::cout << "gfs with " << dg1_gfs_.size() << " dofs generated "<< std::endl; - std::cout << "cc with " << dg1_gfs__cc.size() << " dofs generated "<< std::endl; - - // Set up solution vectors... - using V_R = Dune::PDELab::Backend::Vector<DG1_GFS,DF>; - V_R x_r(dg1_gfs_); - x_r = 0.0; - auto lambda_0000 = [&](const auto& is, const auto& xl){ auto x=is.geometry().global(xl); return exp((-1.0) * ((0.5 - x[2]) * (0.5 - x[2]) + (0.5 - x[1]) * (0.5 - x[1]) + (0.5 - x[0]) * (0.5 - x[0])));; }; - auto func_0000 = Dune::PDELab::makeGridFunctionFromCallable(gv, lambda_0000); - - // Set up (non)linear solvers... - using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_SuperLU; - using SLP = Dune::PDELab::StationaryLinearProblemSolver<GO_r, LinearSolver, V_R>; - LinearSolver ls(false); - double reduction = initree.get<double>("reduction", 1e-12); - SLP slp(go_r, ls, x_r, reduction); - // slp.apply(); - - // Do visualization... - using VTKWriter = Dune::SubsamplingVTKWriter<GV>; - Dune::RefinementIntervals subint(initree.get<int>("vtk.subsamplinglevel", 1)); - VTKWriter vtkwriter(gv, subint); - std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output"); - CuttingPredicate predicate; - Dune::PDELab::addSolutionToVTKWriter(vtkwriter, dg1_gfs_, x_r, Dune::PDELab::vtk::defaultNameScheme(), predicate); - vtkwriter.write(vtkfile, Dune::VTK::ascii); - - // Maybe print residuals and matrices to stdout... - if (initree.get<bool>("printresidual", false)) { - using Dune::PDELab::Backend::native; - V_R r(x_r); - // Interpolate input - auto interpolate_lambda = [] (const auto& x){ - return std::exp(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]); - }; - auto interpolate = Dune::PDELab::makeGridFunctionFromCallable(gv, interpolate_lambda); - Dune::PDELab::interpolate(interpolate,dg1_gfs_,x_r); - Dune::printvector(std::cout, native(x_r), "result of interpolate", "row"); - - r=0.0; - go_r.residual(x_r, r); - Dune::printvector(std::cout, native(r), "residual vector", "row"); - - - std::cout.precision(17); - std::vector<RangeType> residual(16); - for (std::size_t i=0; i<16; ++i){ - residual[i] = native(r)[i]; - } - std::cout << "palpo residual:" << std::endl; - for (std::size_t i=0; i<16; ++i){ - std::cout << residual[i] << " "; - } - std::cout << std::endl; - - - using DGF = Dune::PDELab::DiscreteGridFunction<DG1_GFS,V_R>; - DGF dgf(dg1_gfs_,r); - for (const auto& e : elements(gv)) - { - auto geo = e.geometry(); - DGF::Traits::RangeType eval; - for (int i=0; i<geo.corners(); i++){ - auto local_corner = geo.local(geo.corner(i)); - dgf.evaluate(e, local_corner, eval); - std::cout << "palpo eval: " << eval << std::endl; - } - } - - - - } - if (initree.get<bool>("printmatrix", false)) { - using Dune::PDELab::Backend::native; - V_R r(x_r); - // Interpolate input - auto interpolate_lambda = [] (const auto& x){ - return std::exp(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]); - }; - auto interpolate = Dune::PDELab::makeGridFunctionFromCallable(gv, interpolate_lambda); - Dune::PDELab::interpolate(interpolate,dg1_gfs_,x_r); - using M = typename GO_r::Traits::Jacobian; - M m(go_r); - go_r.jacobian(x_r,m); - using Dune::PDELab::Backend::native; - Dune::printmatrix(std::cout, native(m),"global stiffness matrix","row",9,1); - } - - // // Maybe calculate errors for test results... - // using DG1_GFS__DGF = Dune::PDELab::DiscreteGridFunction<decltype(dg1_gfs_),decltype(x_r)>; - // DG1_GFS__DGF dg1_gfs__dgf(dg1_gfs_,x_r); - // using DifferenceSquaredAdapter_ = Dune::PDELab::DifferenceSquaredAdapter<decltype(func_0000), decltype(dg1_gfs__dgf)>; - // DifferenceSquaredAdapter_ dsa_(func_0000, dg1_gfs__dgf); - // Dune::FieldVector<RangeType, 1> l2error(0.0); - // { - // // L2 error squared of difference between numerical - // // solution and the interpolation of exact solution - // // for treepath () - // typename decltype(dsa_)::Traits::RangeType err(0.0); - // Dune::PDELab::integrateGridFunction(dsa_, err, 10); - - // l2error += err; - // if (gv.comm().rank() == 0){ - // std::cout << "L2 Error for treepath : " << err << std::endl; - // }} - // bool testfail(false); - // using std::abs; - // using std::isnan; - // if (gv.comm().rank() == 0){ - // std::cout << "\nl2errorsquared: " << l2error << std::endl << std::endl; - // } - // if (isnan(l2error[0]) or abs(l2error[0])>1e-4) - // testfail = true; - // return testfail; - return 1; - - } - catch (Dune::Exception& e) - { std::cerr << "Dune reported error: " << e << std::endl; - return 1; - } - catch (std::exception& e) - { std::cerr << "Unknown exception thrown!" << std::endl; - return 1; - } -} diff --git a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_gmsh_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_gmsh_3d_facedir_facemod_variation_driver.cc index 27abe51fefa27120a4c6423639baa7c25a1da543..311c39cdafd074d7fa7e1db464e9f06417d22738 100644 --- a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_gmsh_3d_facedir_facemod_variation_driver.cc +++ b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_gmsh_3d_facedir_facemod_variation_driver.cc @@ -46,12 +46,12 @@ void test_grid(GV gv){ // iterate over all entities of the grid for (const auto& e : elements(gv)) { - std::cout << "## palpo New Element!" << std::endl; + std::cout << "## New Element!" << std::endl; auto geo = e.geometry(); for (int i=0; i<geo.corners(); i++){ auto global_corner = geo.corner(i); auto local_corner = geo.local(global_corner); - std::cout << "palpo global_corner: "; + std::cout << "global_corner: "; for (int i=0; i<3; ++i) std::cout << global_corner[i] << " "; std::cout << std::endl; @@ -59,7 +59,7 @@ void test_grid(GV gv){ for (int i=0; i<geo.corners(); i++){ auto global_corner = geo.corner(i); auto local_corner = geo.local(global_corner); - std::cout << "palpo local_corner: "; + std::cout << "local_corner: "; for (int i=0; i<3; ++i) std::cout << local_corner[i] << " "; std::cout << std::endl; @@ -211,7 +211,7 @@ int main(int argc, char** argv){ native(x_s)[i] = native(x_n)[i]; // std::cout.precision(17); - std::cout << "palpo coefficients:" << std::endl; + std::cout << "Coefficients:" << std::endl; for (std::size_t i=0; i<16; ++i){ std::cout << native(x_s)[i] << ", "; } @@ -228,30 +228,11 @@ int main(int argc, char** argv){ for (std::size_t i=0; i<16; ++i){ residual[i] = native(r)[i]; } - std::cout << "palpo residual:" << std::endl; + std::cout << "Residual:" << std::endl; for (std::size_t i=0; i<16; ++i){ std::cout << residual[i] << " "; } std::cout << std::endl; - - std::vector<RangeType> solution_1 {-0.003747426404075297, -0.21864656102913654, -0.0032054511665800555, -0.28437017759171157, -0.0018737132020376381, -0.2429341647618537, -0.0016027255832900197, -0.29660333834120023, 0.22916620659136877, 0.0037474264040752831, 0.27616815788636451, 0.0032054511665800512, 0.24632027434093215, 0.0018737132020376344, 0.2908996029052367, 0.0016027255832900187}; - - std::vector<RangeType> solution_2 {-0.22216923038681022, -0.003662938655271696, -0.23743729880973413, -0.0030364756689728193, -0.29395296456019571, -0.0018314693276358439, -0.29973173153779054, -0.0015182378344864129, 0.0036629386552717003, 0.2309838664360746, 0.0030364756689728349, 0.24202149140895873, 0.001831469327635845, 0.28737146661901281, 0.001518237834486417, 0.29291440083048437}; - - std::vector<RangeType> solution_3 {-0.21436796383145873, -0.23383260730195662, -0.28383755411465533, -0.29526590985481704, -0.0045507639649666034, -0.0034803883238202818, -0.0036071199470257055, -0.0024060631441813391, 0.22702984288940109, 0.0045507639649666042, 0.27429663847448937, 0.0036071199470257059, 0.23989650875067281, 0.0034803883238202818, 0.28608104498832443, 0.0024060631441813391}; - - std::vector<RangeType> solution_4 {2.3589302088734218e-17, 6.3207344446901369e-18, 6.32073444469013e-18, 1.693635690026303e-18, -0.24497506595070379, -0.29846897291910235, -0.25528824297544883, -0.30346218805770897, 0.24497506595070367, -3.8844456896567597e-17, 0.2984689729191024, -1.0408340855860846e-17, 0.25528824297544883, -1.0408340855860843e-17, 0.30346218805770903, -2.7889065268757736e-18}; - - - std::cout << "palpo is_permuation to solution_1: " - << fuzzy_is_permutation(solution_1, residual) << std::endl; - std::cout << "palpo is_permuation to solution_2: " - << fuzzy_is_permutation(solution_2, residual) << std::endl; - std::cout << "palpo is_permuation to solution_3: " - << fuzzy_is_permutation(solution_3, residual) << std::endl; - std::cout << "palpo is_permuation to solution_4: " - << fuzzy_is_permutation(solution_4, residual) << std::endl; - } if (initree.get<bool>("printmatrix", false)) { using Dune::PDELab::Backend::native; @@ -295,7 +276,7 @@ int main(int argc, char** argv){ // if (isnan(l2error[0]) or abs(l2error[0])>1e-4) // testfail = true; // return testfail; - return 1; + return 0; } catch (Dune::Exception& e) diff --git a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_structured_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_structured_3d_facedir_facemod_variation_driver.cc index 6ddd7fff3613bd6e098bef1065608490179d98b2..53c18b31cedc62cce65b4f0ad7cf74e32f16bcf7 100644 --- a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_structured_3d_facedir_facemod_variation_driver.cc +++ b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_structured_3d_facedir_facemod_variation_driver.cc @@ -113,28 +113,11 @@ int main(int argc, char** argv){ for (std::size_t i=0; i<16; ++i){ residual[i] = native(r)[i]; } - std::cout << "palpo residual:" << std::endl; + std::cout << "residual:" << std::endl; for (std::size_t i=0; i<16; ++i){ std::cout << residual[i] << " "; } std::cout << std::endl; - - - - using DGF = Dune::PDELab::DiscreteGridFunction<DG1_GFS,V_R>; - DGF dgf(dg1_gfs_,r); - for (const auto& e : elements(gv)) - { - auto geo = e.geometry(); - DGF::Traits::RangeType eval; - for (int i=0; i<geo.corners(); i++){ - auto local_corner = geo.local(geo.corner(i)); - dgf.evaluate(e, local_corner, eval); - std::cout << "palpo eval: " << eval << std::endl; - } - } - - } if (initree.get<bool>("printmatrix", false)) { using Dune::PDELab::Backend::native; diff --git a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_unstructured_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_unstructured_3d_facedir_facemod_variation_driver.cc index e68ecd5d2d852ea52c6771ddce91f6b860326524..7c81d42c99b7b965bc237386fe73adc24d50531a 100644 --- a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_unstructured_3d_facedir_facemod_variation_driver.cc +++ b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_unstructured_3d_facedir_facemod_variation_driver.cc @@ -26,16 +26,16 @@ void test_gmsh_grid(GV gv){ // iterate over all entities of the grid for (const auto& e : elements(gv)) { - std::cout << "## palpo New Element!" << std::endl; + std::cout << "## New Element!" << std::endl; auto geo = e.geometry(); for (int i=0; i<geo.corners(); i++){ auto global_corner = geo.corner(i); auto local_corner = geo.local(global_corner); - std::cout << "palpo global_corner: "; + std::cout << "global_corner: "; for (int i=0; i<3; ++i) std::cout << global_corner[i] << " "; std::cout << std::endl; - std::cout << "palpo local_corner: "; + std::cout << "local_corner: "; for (int i=0; i<3; ++i) std::cout << local_corner[i] << " "; std::cout << std::endl; @@ -62,9 +62,10 @@ int main(int argc, char** argv){ GV gv = grid->leafGridView(); test_gmsh_grid(gv); - // Write a gmsh file. Should be the same as grid_06.msh - Dune::GmshWriter<GV> gmsh_writer(gv,6); - gmsh_writer.write("palpo_grid.msh"); + + // // Write a gmsh file. Should be the same as grid_06.msh + // Dune::GmshWriter<GV> gmsh_writer(gv,6); + // gmsh_writer.write("exported_grid.msh"); // Set up finite element maps... using DG1_FEM = Dune::PDELab::QkDGLocalFiniteElementMap<DF, RangeType, 1, 3>; diff --git a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc index f775f507a1632633786fc6f2c0e6ef5c0b9cbb73..520d6e455e86ef62c17db33fd6bdb9b71706b0da 100644 --- a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc +++ b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation_driver.cc @@ -117,12 +117,12 @@ void test_grid(GV gv){ // iterate over all entities of the grid for (const auto& e : elements(gv)) { - std::cout << "## palpo New Element!" << std::endl; + std::cout << "## New Element!" << std::endl; auto geo = e.geometry(); for (int i=0; i<geo.corners(); i++){ auto global_corner = geo.corner(i); auto local_corner = geo.local(global_corner); - std::cout << "palpo global_corner: "; + std::cout << "global_corner: "; for (int i=0; i<3; ++i) std::cout << global_corner[i] << " "; std::cout << std::endl; @@ -130,7 +130,7 @@ void test_grid(GV gv){ for (int i=0; i<geo.corners(); i++){ auto global_corner = geo.corner(i); auto local_corner = geo.local(global_corner); - std::cout << "palpo local_corner: "; + std::cout << "local_corner: "; for (int i=0; i<3; ++i) std::cout << local_corner[i] << " "; std::cout << std::endl; @@ -246,7 +246,7 @@ int main(int argc, char** argv){ for (std::size_t i=0; i<16; ++i){ residual[i] = native(r)[i]; } - std::cout << "palpo residual:" << std::endl; + std::cout << "residual:" << std::endl; for (std::size_t i=0; i<16; ++i){ std::cout << residual[i] << ", "; } @@ -255,7 +255,7 @@ int main(int argc, char** argv){ // One 'correct' numerical solution std::vector<RangeType> solution {-0.0045507639649666424, -0.21358362254819624, -0.0036071199470256972, -0.28388435136052759, -0.00348038832382031, -0.23387940454782871, -0.0024060631441813422, -0.29616724425275981, 0.22624550160613868, 0.0045507639649666424, 0.27434343572036163, 0.0036071199470256972, 0.23994330599654493, 0.00348038832382031, 0.28698237938626719, 0.0024060631441813422}; - std::cout << "palpo is_permuation to solution_1: " + std::cout << "is_permuation: " << fuzzy_is_permutation(solution, residual) << std::endl; if (!fuzzy_is_permutation(solution, residual)){ return 1; diff --git a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_structured_3d_facedir_facemod_variation_driver.cc b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_structured_3d_facedir_facemod_variation_driver.cc index f2eac8227283b73298500f76f3752f63cca7a5e2..7f7f42927ad1fe8c6d4dad0d723f11419b13e8fb 100644 --- a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_structured_3d_facedir_facemod_variation_driver.cc +++ b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_structured_3d_facedir_facemod_variation_driver.cc @@ -136,7 +136,7 @@ int main(int argc, char** argv){ residual[i] = native(r)[i]; } - std::cout << "palpo is_permuation: " + std::cout << "Result of is_permuation: " << fuzzy_is_permutation(solution, residual) << std::endl; if (!fuzzy_is_permutation(solution, residual)) return 1;