/home/docs/checkouts/readthedocs.org/user_builds/ginkgo-test/checkouts/latest/build/doc/doxygen/examples/external-lib-interfacing.hpp Source File

/home/docs/checkouts/readthedocs.org/user_builds/ginkgo-test/checkouts/latest/build/doc/doxygen/examples/external-lib-interfacing.hpp Source File#

Reference API: /home/docs/checkouts/readthedocs.org/user_builds/ginkgo-test/checkouts/latest/build/doc/doxygen/examples/external-lib-interfacing.hpp Source File
Reference API
external-lib-interfacing.hpp
1
1840 *
1841 * /* ---------------------------------------------------------------------
1842 * *
1843 * * This file has been taken verbatim from the deal.ii (version 9.0)
1844 * * examples directory and modified.
1845 * *
1846 * * This example aims to demonstrate the ease with which Ginkgo can
1847 * * be interfaced with other libraries. The only modification/ addition
1848 * * has been to the AdvectionProblem::solve () function.
1849 * *
1850 * */
1851 *
1852 * #include <deal.II/base/function.h>
1853 * #include <deal.II/base/logstream.h>
1854 * #include <deal.II/base/quadrature_lib.h>
1855 * #include <deal.II/dofs/dof_accessor.h>
1856 * #include <deal.II/dofs/dof_handler.h>
1857 * #include <deal.II/dofs/dof_tools.h>
1858 * #include <deal.II/fe/fe_q.h>
1859 * #include <deal.II/fe/fe_values.h>
1860 * #include <deal.II/grid/grid_generator.h>
1861 * #include <deal.II/grid/grid_out.h>
1862 * #include <deal.II/grid/grid_refinement.h>
1863 * #include <deal.II/grid/tria.h>
1864 * #include <deal.II/grid/tria_accessor.h>
1865 * #include <deal.II/grid/tria_iterator.h>
1866 * #include <deal.II/lac/constraint_matrix.h>
1867 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
1868 * #include <deal.II/lac/full_matrix.h>
1869 * #include <deal.II/lac/precondition.h>
1870 * #include <deal.II/lac/solver_bicgstab.h>
1871 * #include <deal.II/lac/sparse_matrix.h>
1872 * #include <deal.II/lac/vector.h>
1873 * #include <deal.II/numerics/data_out.h>
1874 * #include <deal.II/numerics/matrix_tools.h>
1875 * #include <deal.II/numerics/vector_tools.h>
1876 *
1877 * #include <deal.II/base/multithread_info.h>
1878 * #include <deal.II/base/work_stream.h>
1879 *
1880 * #include <deal.II/base/tensor_function.h>
1881 * #include <deal.II/numerics/error_estimator.h>
1882 *
1883 * #include <ginkgo/ginkgo.hpp>
1884 *
1885 * #include <fstream>
1886 * #include <iostream>
1887 *
1888 *
1889 * namespace Step9 {
1890 * using namespace dealii;
1891 *
1892 *
1893 * template <int dim>
1894 * class AdvectionProblem {
1895 * public:
1896 * AdvectionProblem();
1897 * ~AdvectionProblem();
1898 * void run();
1899 *
1900 * private:
1901 * void setup_system();
1902 *
1903 * struct AssemblyScratchData {
1904 * AssemblyScratchData(const FiniteElement<dim>& fe);
1905 * AssemblyScratchData(const AssemblyScratchData& scratch_data);
1906 *
1907 * FEValues<dim> fe_values;
1908 * FEFaceValues<dim> fe_face_values;
1909 * };
1910 *
1911 * struct AssemblyCopyData {
1912 * FullMatrix<double> cell_matrix;
1913 * Vector<double> cell_rhs;
1914 * std::vector<types::global_dof_index> local_dof_indices;
1915 * };
1916 *
1917 * void assemble_system();
1918 * void local_assemble_system(
1919 * const typename DoFHandler<dim>::active_cell_iterator& cell,
1920 * AssemblyScratchData& scratch, AssemblyCopyData& copy_data);
1921 * void copy_local_to_global(const AssemblyCopyData& copy_data);
1922 *
1923 *
1924 * void solve();
1925 * void refine_grid();
1926 * void output_results(const unsigned int cycle) const;
1927 *
1928 * Triangulation<dim> triangulation;
1929 * DoFHandler<dim> dof_handler;
1930 *
1931 * FE_Q<dim> fe;
1932 *
1933 * ConstraintMatrix hanging_node_constraints;
1934 *
1935 * SparsityPattern sparsity_pattern;
1936 * SparseMatrix<double> system_matrix;
1937 *
1938 * Vector<double> solution;
1939 * Vector<double> system_rhs;
1940 * };
1941 *
1942 *
1943 *
1944 * template <int dim>
1945 * class AdvectionField : public TensorFunction<1, dim> {
1946 * public:
1947 * AdvectionField() : TensorFunction<1, dim>() {}
1948 *
1949 * virtual Tensor<1, dim> value(const Point<dim>& p) const;
1950 *
1951 * virtual void value_list(const std::vector<Point<dim>>& points,
1952 * std::vector<Tensor<1, dim>>& values) const;
1953 *
1954 * DeclException2(ExcDimensionMismatch, unsigned int, unsigned int,
1955 * << "The vector has size " << arg1 << " but should have "
1956 * << arg2 << " elements.");
1957 * };
1958 *
1959 *
1960 * template <int dim>
1961 * Tensor<1, dim> AdvectionField<dim>::value(const Point<dim>& p) const
1962 * {
1963 * Point<dim> value;
1964 * value[0] = 2;
1965 * for (unsigned int i = 1; i < dim; ++i)
1966 * value[i] = 1 + 0.8 * std::sin(8 * numbers::PI * p[0]);
1967 *
1968 * return value;
1969 * }
1970 *
1971 *
1972 * template <int dim>
1973 * void AdvectionField<dim>::value_list(const std::vector<Point<dim>>& points,
1974 * std::vector<Tensor<1, dim>>& values) const
1975 * {
1976 * Assert(values.size() == points.size(),
1977 * ExcDimensionMismatch(values.size(), points.size()));
1978 *
1979 * for (unsigned int i = 0; i < points.size(); ++i)
1980 * values[i] = AdvectionField<dim>::value(points[i]);
1981 * }
1982 *
1983 *
1984 * template <int dim>
1985 * class RightHandSide : public Function<dim> {
1986 * public:
1987 * RightHandSide() : Function<dim>() {}
1988 *
1989 * virtual double value(const Point<dim>& p,
1990 * const unsigned int component = 0) const;
1991 *
1992 * virtual void value_list(const std::vector<Point<dim>>& points,
1993 * std::vector<double>& values,
1994 * const unsigned int component = 0) const;
1995 *
1996 * private:
1997 * static const Point<dim> center_point;
1998 * };
1999 *
2000 *
2001 * template <>
2002 * const Point<1> RightHandSide<1>::center_point = Point<1>(-0.75);
2003 *
2004 * template <>
2005 * const Point<2> RightHandSide<2>::center_point = Point<2>(-0.75, -0.75);
2006 *
2007 * template <>
2008 * const Point<3> RightHandSide<3>::center_point = Point<3>(-0.75, -0.75, -0.75);
2009 *
2010 *
2011 * template <int dim>
2012 * double RightHandSide<dim>::value(const Point<dim>& p,
2013 * const unsigned int component) const
2014 * {
2015 * (void)component;
2016 * Assert(component == 0, ExcIndexRange(component, 0, 1));
2017 * const double diameter = 0.1;
2018 * return ((p - center_point).norm_square() < diameter * diameter
2019 * ? .1 / std::pow(diameter, dim)
2020 * : 0);
2021 * }
2022 *
2023 *
2024 * template <int dim>
2025 * void RightHandSide<dim>::value_list(const std::vector<Point<dim>>& points,
2026 * std::vector<double>& values,
2027 * const unsigned int component) const
2028 * {
2029 * Assert(values.size() == points.size(),
2030 * ExcDimensionMismatch(values.size(), points.size()));
2031 *
2032 * for (unsigned int i = 0; i < points.size(); ++i)
2033 * values[i] = RightHandSide<dim>::value(points[i], component);
2034 * }
2035 *
2036 *
2037 * template <int dim>
2038 * class BoundaryValues : public Function<dim> {
2039 * public:
2040 * BoundaryValues() : Function<dim>() {}
2041 *
2042 * virtual double value(const Point<dim>& p,
2043 * const unsigned int component = 0) const;
2044 *
2045 * virtual void value_list(const std::vector<Point<dim>>& points,
2046 * std::vector<double>& values,
2047 * const unsigned int component = 0) const;
2048 * };
2049 *
2050 *
2051 * template <int dim>
2052 * double BoundaryValues<dim>::value(const Point<dim>& p,
2053 * const unsigned int component) const
2054 * {
2055 * (void)component;
2056 * Assert(component == 0, ExcIndexRange(component, 0, 1));
2057 *
2058 * const double sine_term =
2059 * std::sin(16 * numbers::PI * std::sqrt(p.norm_square()));
2060 * const double weight = std::exp(-5 * p.norm_square()) / std::exp(-5.);
2061 * return sine_term * weight;
2062 * }
2063 *
2064 *
2065 * template <int dim>
2066 * void BoundaryValues<dim>::value_list(const std::vector<Point<dim>>& points,
2067 * std::vector<double>& values,
2068 * const unsigned int component) const
2069 * {
2070 * Assert(values.size() == points.size(),
2071 * ExcDimensionMismatch(values.size(), points.size()));
2072 *
2073 * for (unsigned int i = 0; i < points.size(); ++i)
2074 * values[i] = BoundaryValues<dim>::value(points[i], component);
2075 * }
2076 *
2077 *
2078 *
2079 * class GradientEstimation {
2080 * public:
2081 * template <int dim>
2082 * static void estimate(const DoFHandler<dim>& dof,
2083 * const Vector<double>& solution,
2084 * Vector<float>& error_per_cell);
2085 *
2086 * DeclException2(ExcInvalidVectorLength, int, int,
2087 * << "Vector has length " << arg1 << ", but should have "
2088 * << arg2);
2089 * DeclException0(ExcInsufficientDirections);
2090 *
2091 * private:
2092 * template <int dim>
2093 * struct EstimateScratchData {
2094 * EstimateScratchData(const FiniteElement<dim>& fe,
2095 * const Vector<double>& solution,
2096 * Vector<float>& error_per_cell);
2097 * EstimateScratchData(const EstimateScratchData& data);
2098 *
2099 * FEValues<dim> fe_midpoint_value;
2100 * const Vector<double>& solution;
2101 * Vector<float>& error_per_cell;
2102 * };
2103 *
2104 * struct EstimateCopyData {};
2105 *
2106 * template <int dim>
2107 * static void estimate_cell(
2108 * const typename DoFHandler<dim>::active_cell_iterator& cell,
2109 * EstimateScratchData<dim>& scratch_data,
2110 * const EstimateCopyData& copy_data);
2111 * };
2112 *
2113 *
2114 *
2115 *
2116 * template <int dim>
2117 * AdvectionProblem<dim>::AdvectionProblem() : dof_handler(triangulation), fe(1)
2118 * {}
2119 *
2120 *
2121 * template <int dim>
2122 * AdvectionProblem<dim>::~AdvectionProblem()
2123 * {
2124 * dof_handler.clear();
2125 * }
2126 *
2127 *
2128 * template <int dim>
2129 * void AdvectionProblem<dim>::setup_system()
2130 * {
2131 * dof_handler.distribute_dofs(fe);
2132 * hanging_node_constraints.clear();
2133 * DoFTools::make_hanging_node_constraints(dof_handler,
2134 * hanging_node_constraints);
2135 * hanging_node_constraints.close();
2136 *
2137 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
2138 * DoFTools::make_sparsity_pattern(dof_handler, dsp, hanging_node_constraints,
2139 * /*keep_constrained_dofs = */ true);
2140 * sparsity_pattern.copy_from(dsp);
2141 *
2142 * system_matrix.reinit(sparsity_pattern);
2143 *
2144 * solution.reinit(dof_handler.n_dofs());
2145 * system_rhs.reinit(dof_handler.n_dofs());
2146 * }
2147 *
2148 *
2149 * template <int dim>
2150 * void AdvectionProblem<dim>::assemble_system()
2151 * {
2152 * WorkStream::run(dof_handler.begin_active(), dof_handler.end(), *this,
2153 * &AdvectionProblem::local_assemble_system,
2154 * &AdvectionProblem::copy_local_to_global,
2155 * AssemblyScratchData(fe), AssemblyCopyData());
2156 *
2157 *
2158 * hanging_node_constraints.condense(system_matrix);
2159 * hanging_node_constraints.condense(system_rhs);
2160 * }
2161 *
2162 *
2163 * template <int dim>
2164 * AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
2165 * const FiniteElement<dim>& fe)
2166 * : fe_values(fe, QGauss<dim>(2),
2167 * update_values | update_gradients | update_quadrature_points |
2168 * update_JxW_values),
2169 * fe_face_values(fe, QGauss<dim - 1>(2),
2170 * update_values | update_quadrature_points |
2171 * update_JxW_values | update_normal_vectors)
2172 * {}
2173 *
2174 *
2175 * template <int dim>
2176 * AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
2177 * const AssemblyScratchData& scratch_data)
2178 * : fe_values(scratch_data.fe_values.get_fe(),
2179 * scratch_data.fe_values.get_quadrature(),
2180 * update_values | update_gradients | update_quadrature_points |
2181 * update_JxW_values),
2182 * fe_face_values(scratch_data.fe_face_values.get_fe(),
2183 * scratch_data.fe_face_values.get_quadrature(),
2184 * update_values | update_quadrature_points |
2185 * update_JxW_values | update_normal_vectors)
2186 * {}
2187 *
2188 *
2189 * template <int dim>
2190 * void AdvectionProblem<dim>::local_assemble_system(
2191 * const typename DoFHandler<dim>::active_cell_iterator& cell,
2192 * AssemblyScratchData& scratch_data, AssemblyCopyData& copy_data)
2193 * {
2194 * const AdvectionField<dim> advection_field;
2195 * const RightHandSide<dim> right_hand_side;
2196 * const BoundaryValues<dim> boundary_values;
2197 *
2198 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
2199 * const unsigned int n_q_points =
2200 * scratch_data.fe_values.get_quadrature().size();
2201 * const unsigned int n_face_q_points =
2202 * scratch_data.fe_face_values.get_quadrature().size();
2203 *
2204 * copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
2205 * copy_data.cell_rhs.reinit(dofs_per_cell);
2206 *
2207 * copy_data.local_dof_indices.resize(dofs_per_cell);
2208 *
2209 * std::vector<double> rhs_values(n_q_points);
2210 * std::vector<Tensor<1, dim>> advection_directions(n_q_points);
2211 * std::vector<double> face_boundary_values(n_face_q_points);
2212 * std::vector<Tensor<1, dim>> face_advection_directions(n_face_q_points);
2213 *
2214 *
2215 * scratch_data.fe_values.reinit(cell);
2216 *
2217 * advection_field.value_list(scratch_data.fe_values.get_quadrature_points(),
2218 * advection_directions);
2219 * right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(),
2220 * rhs_values);
2221 *
2222 * const double delta = 0.1 * cell->diameter();
2223 *
2224 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2225 * for (unsigned int i = 0; i < dofs_per_cell; ++i) {
2226 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2227 * copy_data.cell_matrix(i, j) +=
2228 * ((advection_directions[q_point] *
2229 * scratch_data.fe_values.shape_grad(j, q_point) *
2230 * (scratch_data.fe_values.shape_value(i, q_point) +
2231 * delta *
2232 * (advection_directions[q_point] *
2233 * scratch_data.fe_values.shape_grad(i, q_point)))) *
2234 * scratch_data.fe_values.JxW(q_point));
2235 *
2236 * copy_data.cell_rhs(i) +=
2237 * ((scratch_data.fe_values.shape_value(i, q_point) +
2238 * delta * (advection_directions[q_point] *
2239 * scratch_data.fe_values.shape_grad(i, q_point))) *
2240 * rhs_values[q_point] * scratch_data.fe_values.JxW(q_point));
2241 * }
2242 *
2243 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2244 * ++face)
2245 * if (cell->face(face)->at_boundary()) {
2246 * scratch_data.fe_face_values.reinit(cell, face);
2247 *
2248 * boundary_values.value_list(
2249 * scratch_data.fe_face_values.get_quadrature_points(),
2250 * face_boundary_values);
2251 * advection_field.value_list(
2252 * scratch_data.fe_face_values.get_quadrature_points(),
2253 * face_advection_directions);
2254 *
2255 * for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
2256 * if (scratch_data.fe_face_values.normal_vector(q_point) *
2257 * face_advection_directions[q_point] <
2258 * 0)
2259 * for (unsigned int i = 0; i < dofs_per_cell; ++i) {
2260 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2261 * copy_data.cell_matrix(i, j) -=
2262 * (face_advection_directions[q_point] *
2263 * scratch_data.fe_face_values.normal_vector(
2264 * q_point) *
2265 * scratch_data.fe_face_values.shape_value(
2266 * i, q_point) *
2267 * scratch_data.fe_face_values.shape_value(
2268 * j, q_point) *
2269 * scratch_data.fe_face_values.JxW(q_point));
2270 *
2271 * copy_data.cell_rhs(i) -=
2272 * (face_advection_directions[q_point] *
2273 * scratch_data.fe_face_values.normal_vector(
2274 * q_point) *
2275 * face_boundary_values[q_point] *
2276 * scratch_data.fe_face_values.shape_value(i,
2277 * q_point) *
2278 * scratch_data.fe_face_values.JxW(q_point));
2279 * }
2280 * }
2281 *
2282 *
2283 * cell->get_dof_indices(copy_data.local_dof_indices);
2284 * }
2285 *
2286 *
2287 * template <int dim>
2288 * void AdvectionProblem<dim>::copy_local_to_global(
2289 * const AssemblyCopyData& copy_data)
2290 * {
2291 * for (unsigned int i = 0; i < copy_data.local_dof_indices.size(); ++i) {
2292 * for (unsigned int j = 0; j < copy_data.local_dof_indices.size(); ++j)
2293 * system_matrix.add(copy_data.local_dof_indices[i],
2294 * copy_data.local_dof_indices[j],
2295 * copy_data.cell_matrix(i, j));
2296 *
2297 * system_rhs(copy_data.local_dof_indices[i]) += copy_data.cell_rhs(i);
2298 * }
2299 * }
2300 *
2301 * template <int dim>
2302 * void AdvectionProblem<dim>::solve()
2303 * {
2304 * Assert(system_matrix.m() == system_matrix.n(), ExcNotQuadratic());
2305 * auto num_rows = system_matrix.m();
2306 *
2307 * std::vector<double> rhs(num_rows);
2308 * std::copy(system_rhs.begin(), system_rhs.begin() + num_rows, rhs.begin());
2309 *
2310 * using vec = gko::matrix::Dense<>;
2311 * using mtx = gko::matrix::Csr<>;
2312 * using bicgstab = gko::solver::Bicgstab<>;
2313 * using bj = gko::preconditioner::Jacobi<>;
2314 * using val_array = gko::array<double>;
2315 *
2316 * std::shared_ptr<gko::Executor> exec = gko::ReferenceExecutor::create();
2317 *
2318 * auto b = vec::create(exec, gko::dim<2>(num_rows, 1),
2319 * val_array::view(exec, num_rows, rhs.data()), 1);
2320 * auto x = vec::create(exec, gko::dim<2>(num_rows, 1));
2321 * auto A = mtx::create(exec, gko::dim<2>(num_rows),
2322 * system_matrix.n_nonzero_elements());
2323 * mtx::value_type* values = A->get_values();
2324 * mtx::index_type* row_ptr = A->get_row_ptrs();
2325 * mtx::index_type* col_idx = A->get_col_idxs();
2326 *
2327 * row_ptr[0] = 0;
2328 * for (auto row = 1; row <= num_rows; ++row) {
2329 * row_ptr[row] = row_ptr[row - 1] + system_matrix.get_row_length(row - 1);
2330 * }
2331 *
2332 * std::vector<mtx::index_type> ptrs(num_rows + 1);
2333 * std::copy(A->get_row_ptrs(), A->get_row_ptrs() + num_rows + 1,
2334 * ptrs.begin());
2335 * for (auto row = 0; row < system_matrix.m(); ++row) {
2336 * for (auto p = system_matrix.begin(row); p != system_matrix.end(row);
2337 * ++p) {
2338 * col_idx[ptrs[row]] = p->column();
2339 * values[ptrs[row]] = p->value();
2340 *
2341 * ++ptrs[row];
2342 * }
2343 * }
2344 *
2345 * auto solver_gen =
2346 * bicgstab::build()
2347 * .with_criteria(
2348 * gko::stop::Iteration::build().with_max_iters(1000),
2349 * gko::stop::ResidualNorm<>::build().with_reduction_factor(1e-12))
2350 * .with_preconditioner(bj::build())
2351 * .on(exec);
2352 * auto solver = solver_gen->generate(gko::give(A));
2353 *
2354 * solver->apply(b, x);
2355 *
2356 * std::copy(x->get_values(), x->get_values() + num_rows, solution.begin());
2357 *
2358 * /******************************************************
2359 * * deal.ii internal solver. Here for reference.
2360 * SolverControl solver_control (1000, 1e-12);
2361 * SolverBicgstab<> bicgstab (solver_control);
2362 *
2363 * PreconditionJacobi<> preconditioner;
2364 * preconditioner.initialize(system_matrix, 1.0);
2365 *
2366 * bicgstab.solve (system_matrix, solution, system_rhs,
2367 * preconditioner);
2368 * *******************************************************/
2369 *
2370 * hanging_node_constraints.distribute(solution);
2371 * }
2372 *
2373 *
2374 * template <int dim>
2375 * void AdvectionProblem<dim>::refine_grid()
2376 * {
2377 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2378 *
2379 * GradientEstimation::estimate(dof_handler, solution,
2380 * estimated_error_per_cell);
2381 *
2382 * GridRefinement::refine_and_coarsen_fixed_number(
2383 * triangulation, estimated_error_per_cell, 0.5, 0.03);
2384 *
2385 * triangulation.execute_coarsening_and_refinement();
2386 * }
2387 *
2388 *
2389 * template <int dim>
2390 * void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
2391 * {
2392 * {
2393 * GridOut grid_out;
2394 * std::ofstream output("grid-" + std::to_string(cycle) + ".eps");
2395 * grid_out.write_eps(triangulation, output);
2396 * }
2397 *
2398 * {
2399 * DataOut<dim> data_out;
2400 * data_out.attach_dof_handler(dof_handler);
2401 * data_out.add_data_vector(solution, "solution");
2402 * data_out.build_patches();
2403 *
2404 * std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
2405 * data_out.write_vtk(output);
2406 * }
2407 * }
2408 *
2409 *
2410 * template <int dim>
2411 * void AdvectionProblem<dim>::run()
2412 * {
2413 * for (unsigned int cycle = 0; cycle < 6; ++cycle) {
2414 * std::cout << "Cycle " << cycle << ':' << std::endl;
2415 *
2416 * if (cycle == 0) {
2417 * GridGenerator::hyper_cube(triangulation, -1, 1);
2418 * triangulation.refine_global(4);
2419 * } else {
2420 * refine_grid();
2421 * }
2422 *
2423 *
2424 * std::cout << " Number of active cells: "
2425 * << triangulation.n_active_cells() << std::endl;
2426 *
2427 * setup_system();
2428 *
2429 * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
2430 * << std::endl;
2431 *
2432 * assemble_system();
2433 * solve();
2434 * output_results(cycle);
2435 * }
2436 * }
2437 *
2438 *
2439 *
2440 * template <int dim>
2441 * GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
2442 * const FiniteElement<dim>& fe, const Vector<double>& solution,
2443 * Vector<float>& error_per_cell)
2444 * : fe_midpoint_value(fe, QMidpoint<dim>(),
2445 * update_values | update_quadrature_points),
2446 * solution(solution),
2447 * error_per_cell(error_per_cell)
2448 * {}
2449 *
2450 *
2451 * template <int dim>
2452 * GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
2453 * const EstimateScratchData& scratch_data)
2454 * : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(),
2455 * scratch_data.fe_midpoint_value.get_quadrature(),
2456 * update_values | update_quadrature_points),
2457 * solution(scratch_data.solution),
2458 * error_per_cell(scratch_data.error_per_cell)
2459 * {}
2460 *
2461 *
2462 * template <int dim>
2463 * void GradientEstimation::estimate(const DoFHandler<dim>& dof_handler,
2464 * const Vector<double>& solution,
2465 * Vector<float>& error_per_cell)
2466 * {
2467 * Assert(error_per_cell.size() ==
2468 * dof_handler.get_triangulation().n_active_cells(),
2469 * ExcInvalidVectorLength(
2470 * error_per_cell.size(),
2471 * dof_handler.get_triangulation().n_active_cells()));
2472 *
2473 * WorkStream::run(dof_handler.begin_active(), dof_handler.end(),
2474 * &GradientEstimation::template estimate_cell<dim>,
2475 * std::function<void(const EstimateCopyData&)>(),
2476 * EstimateScratchData<dim>(dof_handler.get_fe(), solution,
2477 * error_per_cell),
2478 * EstimateCopyData());
2479 * }
2480 *
2481 *
2482 * template <int dim>
2483 * void GradientEstimation::estimate_cell(
2484 * const typename DoFHandler<dim>::active_cell_iterator& cell,
2485 * EstimateScratchData<dim>& scratch_data, const EstimateCopyData&)
2486 * {
2487 * Tensor<2, dim> Y;
2488 *
2489 *
2490 * std::vector<typename DoFHandler<dim>::active_cell_iterator>
2491 * active_neighbors;
2492 * active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
2493 * GeometryInfo<dim>::max_children_per_face);
2494 *
2495 * scratch_data.fe_midpoint_value.reinit(cell);
2496 *
2497 * Tensor<1, dim> projected_gradient;
2498 *
2499 *
2500 * active_neighbors.clear();
2501 * for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
2502 * ++face_no)
2503 * if (!cell->at_boundary(face_no)) {
2504 * const typename DoFHandler<dim>::face_iterator face =
2505 * cell->face(face_no);
2506 * const typename DoFHandler<dim>::cell_iterator neighbor =
2507 * cell->neighbor(face_no);
2508 *
2509 * if (neighbor->active())
2510 * active_neighbors.push_back(neighbor);
2511 * else {
2512 * if (dim == 1) {
2513 * typename DoFHandler<dim>::cell_iterator neighbor_child =
2514 * neighbor;
2515 * while (neighbor_child->has_children())
2516 * neighbor_child =
2517 * neighbor_child->child(face_no == 0 ? 1 : 0);
2518 *
2519 * Assert(
2520 * neighbor_child->neighbor(face_no == 0 ? 1 : 0) == cell,
2521 * ExcInternalError());
2522 *
2523 * active_neighbors.push_back(neighbor_child);
2524 * } else
2525 * for (unsigned int subface_no = 0;
2526 * subface_no < face->n_children(); ++subface_no)
2527 * active_neighbors.push_back(
2528 * cell->neighbor_child_on_subface(face_no,
2529 * subface_no));
2530 * }
2531 * }
2532 *
2533 * const Point<dim> this_center =
2534 * scratch_data.fe_midpoint_value.quadrature_point(0);
2535 *
2536 * std::vector<double> this_midpoint_value(1);
2537 * scratch_data.fe_midpoint_value.get_function_values(scratch_data.solution,
2538 * this_midpoint_value);
2539 *
2540 *
2541 * std::vector<double> neighbor_midpoint_value(1);
2542 * typename std::vector<typename DoFHandler<dim>::active_cell_iterator>::
2543 * const_iterator neighbor_ptr = active_neighbors.begin();
2544 * for (; neighbor_ptr != active_neighbors.end(); ++neighbor_ptr) {
2545 * const typename DoFHandler<dim>::active_cell_iterator neighbor =
2546 * *neighbor_ptr;
2547 *
2548 * scratch_data.fe_midpoint_value.reinit(neighbor);
2549 * const Point<dim> neighbor_center =
2550 * scratch_data.fe_midpoint_value.quadrature_point(0);
2551 *
2552 * scratch_data.fe_midpoint_value.get_function_values(
2553 * scratch_data.solution, neighbor_midpoint_value);
2554 *
2555 * Tensor<1, dim> y = neighbor_center - this_center;
2556 * const double distance = y.norm();
2557 * y /= distance;
2558 *
2559 * for (unsigned int i = 0; i < dim; ++i)
2560 * for (unsigned int j = 0; j < dim; ++j) Y[i][j] += y[i] * y[j];
2561 *
2562 * projected_gradient +=
2563 * (neighbor_midpoint_value[0] - this_midpoint_value[0]) / distance *
2564 * y;
2565 * }
2566 *
2567 * AssertThrow(determinant(Y) != 0, ExcInsufficientDirections());
2568 *
2569 * const Tensor<2, dim> Y_inverse = invert(Y);
2570 *
2571 * Tensor<1, dim> gradient = Y_inverse * projected_gradient;
2572 *
2573 * scratch_data.error_per_cell(cell->active_cell_index()) =
2574 * (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) *
2575 * std::sqrt(gradient.norm_square()));
2576 * }
2577 * } // namespace Step9
2578 *
2579 *
2580 *
2581 * int main()
2582 * {
2583 * try {
2584 * dealii::MultithreadInfo::set_thread_limit();
2585 *
2586 * Step9::AdvectionProblem<2> advection_problem_2d;
2587 * advection_problem_2d.run();
2588 * } catch (std::exception& exc) {
2589 * std::cerr << std::endl
2590 * << std::endl
2591 * << "----------------------------------------------------"
2592 * << std::endl;
2593 * std::cerr << "Exception on processing: " << std::endl
2594 * << exc.what() << std::endl
2595 * << "Aborting!" << std::endl
2596 * << "----------------------------------------------------"
2597 * << std::endl;
2598 * return 1;
2599 * } catch (...) {
2600 * std::cerr << std::endl
2601 * << std::endl
2602 * << "----------------------------------------------------"
2603 * << std::endl;
2604 * std::cerr << "Unknown exception!" << std::endl
2605 * << "Aborting!" << std::endl
2606 * << "----------------------------------------------------"
2607 * << std::endl;
2608 * return 1;
2609 * }
2610 *
2611 * return 0;
2612 * }
2613 * @endcode
2614*/
Step9::AdvectionProblem
Definition external-lib-interfacing.cpp:73
gko::array
Definition array.hpp:166
gko::batch::solver::Bicgstab
Definition batch_bicgstab.hpp:50
gko::matrix::Csr
Definition csr.hpp:123
gko::matrix::Dense< ValueType >::create
static std::unique_ptr< Dense > create(std::shared_ptr< const Executor > exec, const dim< 2 > &size={}, size_type stride=0)
gko::preconditioner::Jacobi
Definition jacobi.hpp:189
gko::solver::Bicgstab
Definition bicgstab.hpp:53
gko::stop::ResidualNorm
Definition residual_norm.hpp:113
gko::solver::multigrid::cycle
cycle
Definition multigrid.hpp:54
gko::give
std::remove_reference< OwningPointer >::type && give(OwningPointer &&p)
Definition utils_helper.hpp:247
gko::dim
Definition dim.hpp:26
Generated by doxygen 1.9.8