diff --git a/.noname.gv.swp b/.noname.gv.swp new file mode 100644 index 0000000..7bcefd7 Binary files /dev/null and b/.noname.gv.swp differ diff --git a/libnufi_lib.a b/libnufi_lib.a index da7028d..dd03ba6 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/noname.gv b/noname.gv new file mode 100644 index 0000000..5ea5bf6 --- /dev/null +++ b/noname.gv @@ -0,0 +1,13 @@ +digraph G { + graph [layout=dot rankdir=LR] + +// This is just an example for you to use as a template. +// Edit as you like. Whenever you save a legal graph +// the layout in the graphviz window will be updated. + + vim [href="http://www.vim.org/"] + dot [href="http://www.graphviz.org/"] + vimdot [href="file:///usr/bin/vimdot"] + + {vim dot} -> vimdot +} diff --git a/nufi/cells.h b/nufi/cells.h index e1f74d4..36c61a2 100644 --- a/nufi/cells.h +++ b/nufi/cells.h @@ -19,14 +19,18 @@ template struct CellInfo { double h; }; +template struct CellLocation { + const CellInfo *info; + Point reference_point; +}; + template class CellLocator { public: using CellIterator = typename DoFHandler::active_cell_iterator; void rebuild(const DoFHandler &dof_handler, const Triangulation &triangulation); - - CellIterator locate(const Point &p) const; + CellLocation locate(const Point &p) const; private: std::vector> cells; @@ -58,8 +62,7 @@ void CellLocator::rebuild(const DoFHandler &dof_handler, } template -typename DoFHandler::active_cell_iterator -CellLocator::locate(const Point &p) const { +CellLocation CellLocator::locate(const Point &p) const { static_assert(dim == 1, "Current CellLocator implementation only supports 1D."); @@ -82,7 +85,12 @@ CellLocator::locate(const Point &p) const { AssertThrow(x >= it->lower[0] - 1e-12 && x <= it->upper[0] + 1e-12, ExcMessage("CellLocator failed to find containing cell.")); - return it->cell; + CellLocation location; + + location.info = &(*it); + location.reference_point[0] = (p[0] - it->lower[0]) / it->h; + + return location; } #endif // !CELLS_H diff --git a/nufi/parameters.h b/nufi/parameters.h index bdf74fa..4bcf0be 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -24,8 +24,8 @@ constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV; // deal.ii options constexpr unsigned int GLOBAL_REFINEMENT = 6; -constexpr unsigned int FE_DEGREE = 2; -constexpr unsigned int CONVERGENCE_ITERATIONS = 15000; +constexpr unsigned int FE_DEGREE = 3; +constexpr unsigned int CONVERGENCE_ITERATIONS = 5000; constexpr double CONVERGENCE_LIMIT = 1e-8; constexpr double EPS = 0.01; @@ -35,10 +35,10 @@ constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) // NUFI options constexpr double DT = 1. / 10.; constexpr unsigned int TMAX = 100; -constexpr unsigned int REFINE_FREQUENCY = 5; +constexpr unsigned int REFINE_FREQUENCY = 10; // Plotting options -constexpr int PLOT_FREQUENCY = 5; +constexpr int PLOT_FREQUENCY = 10; constexpr size_t PLOT_NX = CALC_NX; constexpr double PLOT_DX = LX / PLOT_NX; const std::string PLOT_DIR = "results/"; diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index 696f24b..2db18db 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -132,7 +132,7 @@ void PoissonProblem::set_rhs_function( template PoissonProblem::PoissonProblem(unsigned int degree) : triangulation(Triangulation::limit_level_difference_at_vertices), - fe(degree), dof_handler(triangulation), mapping(degree) {} + dof_handler(triangulation), fe(degree), mapping(degree) {} template std::vector @@ -208,12 +208,14 @@ std::vector PoissonProblem::eval_vector_grad( for (unsigned int p = 0; p < points.size(); ++p) { - const auto cell = cell_locator.locate(points[p]); + const auto cell_location = cell_locator.locate(points[p]); - cell->get_dof_values(solution, local_solution_buffer.begin(), - local_solution_buffer.end()); + cell_location.info->cell->get_dof_values( + solution, local_solution_buffer.begin(), local_solution_buffer.end()); - evaluator->reinit(cell, ArrayView>(&points[p], 1)); + evaluator->reinit( + cell_location.info->cell, + ArrayView>(&cell_location.reference_point, 1)); evaluator->evaluate(local_solution_buffer, EvaluationFlags::gradients); values[p] = evaluator->get_gradient(0)[0]; diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 7171552..be078a4 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -206,17 +206,17 @@ void NuFISolver::run() { << std::endl; // START: diagnostics - std::cout << "cells = " << poisson.triangulation.n_active_cells() - << " dofs = " << poisson.dof_handler.n_dofs() << std::endl; - double min_h = 1e100; - double max_h = 0; - - for (auto cell : poisson.triangulation.active_cell_iterators()) { - min_h = std::min(min_h, cell->diameter()); - max_h = std::max(max_h, cell->diameter()); - } - - std::cout << "h ratio = " << max_h / min_h << std::endl; + std::cout << "cells = " << poisson.triangulation.n_active_cells() << "\n" + << " dofs = " << poisson.dof_handler.n_dofs() << "\n"; + // double min_h = 1e100; + // double max_h = 0; + // + // for (auto cell : poisson.triangulation.active_cell_iterators()) { + // min_h = std::min(min_h, cell->diameter()); + // max_h = std::max(max_h, cell->diameter()); + // } + // + // std::cout << "h ratio = " << max_h / min_h << std::endl; // END: diagnostics @@ -248,7 +248,7 @@ void NuFISolver::run() { poisson.coarse_and_refine_grid(it, phi_history); refine_time = timer.elapsed() - refine_start; std::cout << "Refinement step done in " - << std::to_string(std::floor(refine_time)) << "[s]"; + << std::to_string(std::floor(refine_time)) << "[s]" << "\n"; } double timer_elapsed = timer.elapsed();