From: Wolfgang Bangerth Date: Wed, 8 Nov 2017 17:42:04 +0000 (-0700) Subject: Consistently use UNIX line endings. X-Git-Tag: v9.0.0-rc1~802^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a5b92925316db9d254b2e368041e9b6710053cb2;p=dealii.git Consistently use UNIX line endings. --- diff --git a/tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc b/tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc index 9b3b2cb032..e5ae2c5ee8 100644 --- a/tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc +++ b/tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc @@ -1,117 +1,117 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2016 - 2017 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -// Test of basic functionality: -// - Taped doubles -// - Multiple dependent functions and Jacobian computations -// - Update to using C++ data structures where possible -// - Test that the tracing of active dependent/independent variables still -// works with these structures -// Adapted from https://github.com/Homebrew/homebrew-science/blob/master/adol-c.rb - -#include "../tests.h" -#include -#include -#include - -#include - -void test_reset_vector_values (const bool reset_values, const int tape_index) -{ - const unsigned int m = 5; // Dependents - const unsigned int n = 10; // Independents - - std::vector xp (n, 0.0); - std::vector yp (m, 0.0); - std::vector x (n, 1.0); // Dep. variable values initially set here - std::vector y (m, 1.0); - - if (reset_values == false) - for (unsigned int i = 0; i < n; i++) - xp[i] = (i + 1.0) / (2.0 + i); - - trace_on(tape_index); - for (unsigned int i = 0; i < n; i++) - { - x[i] <<= xp[i]; - for (unsigned int j = 0; j < m; ++j) - y[j] *= (j+1)*x[i]; - } - for (unsigned int j = 0; j < m; ++j) - y[j] >>= yp[j]; - - trace_off(); - // tapestats(1, tape_stats); - - // --- Change values --- - if (reset_values == true) - for (unsigned int i = 0; i < n; i++) - xp[i] = (i + 1.0) / (2.0 + i); - - // --- Functions --- - double *f = new double[m]; - function(tape_index, m, n, xp.data(), f); - - deallog << "Evaluation points:" << std::endl; - for (unsigned int i = 0; i < n; ++i) - deallog - << " x[" << i << "]: " << xp[i] - << std::endl; - - deallog << "Function values:" << std::endl; - for (unsigned int j = 0; j < m; ++j) - deallog - << " f[" << j << "]: " << f[j] - << " y[" << j << "]: " << yp[j] - << std::endl; - - // --- Jacobian --- - - double **J = new double*[m]; - for (unsigned int j = 0; j < m; ++j) - J[j] = new double[n]; - - jacobian(tape_index, m, n, xp.data(), J); - - deallog << "Function jacobian J:" << std::endl; - for (unsigned int j = 0; j < m; j++) - { - for (unsigned int i = 0; i < n; i++) - deallog << J[j][i] << (i +#include +#include + +#include + +void test_reset_vector_values (const bool reset_values, const int tape_index) +{ + const unsigned int m = 5; // Dependents + const unsigned int n = 10; // Independents + + std::vector xp (n, 0.0); + std::vector yp (m, 0.0); + std::vector x (n, 1.0); // Dep. variable values initially set here + std::vector y (m, 1.0); + + if (reset_values == false) + for (unsigned int i = 0; i < n; i++) + xp[i] = (i + 1.0) / (2.0 + i); + + trace_on(tape_index); + for (unsigned int i = 0; i < n; i++) + { + x[i] <<= xp[i]; + for (unsigned int j = 0; j < m; ++j) + y[j] *= (j+1)*x[i]; + } + for (unsigned int j = 0; j < m; ++j) + y[j] >>= yp[j]; + + trace_off(); + // tapestats(1, tape_stats); + + // --- Change values --- + if (reset_values == true) + for (unsigned int i = 0; i < n; i++) + xp[i] = (i + 1.0) / (2.0 + i); + + // --- Functions --- + double *f = new double[m]; + function(tape_index, m, n, xp.data(), f); + + deallog << "Evaluation points:" << std::endl; + for (unsigned int i = 0; i < n; ++i) + deallog + << " x[" << i << "]: " << xp[i] + << std::endl; + + deallog << "Function values:" << std::endl; + for (unsigned int j = 0; j < m; ++j) + deallog + << " f[" << j << "]: " << f[j] + << " y[" << j << "]: " << yp[j] + << std::endl; + + // --- Jacobian --- + + double **J = new double*[m]; + for (unsigned int j = 0; j < m; ++j) + J[j] = new double[n]; + + jacobian(tape_index, m, n, xp.data(), J); + + deallog << "Function jacobian J:" << std::endl; + for (unsigned int j = 0; j < m; j++) + { + for (unsigned int i = 0; i < n; i++) + deallog << J[j][i] << (i -#include -#include - -#include - -int main(void) -{ - initlog(); - - const unsigned int m = 5; // Dependents - const unsigned int n = 10; // Independents - std::size_t tape_stats[STAT_SIZE]; - - double *xp = new double[n]; - double *yp = new double[m]; - adouble *x = new adouble[n]; - adouble *y = new adouble[m]; - - for (unsigned int i = 0; i < n; i++) - xp[i] = (i + 1.0) / (2.0 + i); - - for (unsigned int j = 0; j < m; ++j) - y[j] = 1.0; - - trace_on(1); - for (unsigned int i = 0; i < n; i++) - { - x[i] <<= xp[i]; - for (unsigned int j = 0; j < m; ++j) - y[j] *= (j+1)*x[i]; - } - for (unsigned int j = 0; j < m; ++j) - y[j] >>= yp[j]; - - delete[] x; - delete[] y; - - trace_off(); - tapestats(1, tape_stats); - - // --- Functions --- - double *f = new double[m]; - function(1, m, n, xp, f); - - deallog << "Evaluation points:" << std::endl; - for (unsigned int i = 0; i < n; ++i) - deallog - << " x[" << i << "]: " << xp[i] - << std::endl; - - deallog << "Function values:" << std::endl; - for (unsigned int j = 0; j < m; ++j) - deallog - << " f[" << j << "]: " << f[j] - << " y[" << j << "]: " << yp[j] - << std::endl; - - // --- Jacobian --- - - double **J = new double*[m]; - for (unsigned int j = 0; j < m; ++j) - J[j] = new double[n]; - - jacobian(1, m, n, xp, J); - - deallog << "Function jacobian J:" << std::endl; - for (unsigned int j = 0; j < m; j++) - { - for (unsigned int i = 0; i < n; i++) - deallog << J[j][i] << (i +#include +#include + +#include + +int main(void) +{ + initlog(); + + const unsigned int m = 5; // Dependents + const unsigned int n = 10; // Independents + std::size_t tape_stats[STAT_SIZE]; + + double *xp = new double[n]; + double *yp = new double[m]; + adouble *x = new adouble[n]; + adouble *y = new adouble[m]; + + for (unsigned int i = 0; i < n; i++) + xp[i] = (i + 1.0) / (2.0 + i); + + for (unsigned int j = 0; j < m; ++j) + y[j] = 1.0; + + trace_on(1); + for (unsigned int i = 0; i < n; i++) + { + x[i] <<= xp[i]; + for (unsigned int j = 0; j < m; ++j) + y[j] *= (j+1)*x[i]; + } + for (unsigned int j = 0; j < m; ++j) + y[j] >>= yp[j]; + + delete[] x; + delete[] y; + + trace_off(); + tapestats(1, tape_stats); + + // --- Functions --- + double *f = new double[m]; + function(1, m, n, xp, f); + + deallog << "Evaluation points:" << std::endl; + for (unsigned int i = 0; i < n; ++i) + deallog + << " x[" << i << "]: " << xp[i] + << std::endl; + + deallog << "Function values:" << std::endl; + for (unsigned int j = 0; j < m; ++j) + deallog + << " f[" << j << "]: " << f[j] + << " y[" << j << "]: " << yp[j] + << std::endl; + + // --- Jacobian --- + + double **J = new double*[m]; + for (unsigned int j = 0; j < m; ++j) + J[j] = new double[n]; + + jacobian(1, m, n, xp, J); + + deallog << "Function jacobian J:" << std::endl; + for (unsigned int j = 0; j < m; j++) + { + for (unsigned int i = 0; i < n; i++) + deallog << J[j][i] << (i -#include - -#include - -int main(void) -{ - initlog(); - - const unsigned int m = 5; // Dependents - const unsigned int n = 10; // Independents - adtl::setNumDir(n); - - adtl::adouble *x = new adtl::adouble[n]; - for (unsigned int i = 0; i < n; i++) - { - x[i] = (i + 1.0) / (2.0 + i); - x[i].setADValue(i,1); - } - - adtl::adouble *y = new adtl::adouble[m]; - for (unsigned int j = 0; j < m; ++j) - y[j] = 1.0; - - for (unsigned int i = 0; i < n; i++) - for (unsigned int j = 0; j < m; ++j) - y[j] *= (j+1)*x[i]; - - // --- Functions --- - - deallog << "Evaluation points:" << std::endl; - for (unsigned int i = 0; i < n; ++i) - deallog - << " x[" << i << "]: " << x[i].getValue() - << std::endl; - - deallog << "Function values:" << std::endl; - for (unsigned int j = 0; j < m; ++j) - deallog - << " y[" << j << "]: " << y[j].getValue() - << std::endl; - - // --- Jacobian --- - - deallog << "Function jacobian J:" << std::endl; - for (unsigned int j = 0; j < m; j++) - { - for (unsigned int i = 0; i < n; i++) - deallog << y[j].getADValue(i) << (i +#include + +#include + +int main(void) +{ + initlog(); + + const unsigned int m = 5; // Dependents + const unsigned int n = 10; // Independents + adtl::setNumDir(n); + + adtl::adouble *x = new adtl::adouble[n]; + for (unsigned int i = 0; i < n; i++) + { + x[i] = (i + 1.0) / (2.0 + i); + x[i].setADValue(i,1); + } + + adtl::adouble *y = new adtl::adouble[m]; + for (unsigned int j = 0; j < m; ++j) + y[j] = 1.0; + + for (unsigned int i = 0; i < n; i++) + for (unsigned int j = 0; j < m; ++j) + y[j] *= (j+1)*x[i]; + + // --- Functions --- + + deallog << "Evaluation points:" << std::endl; + for (unsigned int i = 0; i < n; ++i) + deallog + << " x[" << i << "]: " << x[i].getValue() + << std::endl; + + deallog << "Function values:" << std::endl; + for (unsigned int j = 0; j < m; ++j) + deallog + << " y[" << j << "]: " << y[j].getValue() + << std::endl; + + // --- Jacobian --- + + deallog << "Function jacobian J:" << std::endl; + for (unsigned int j = 0; j < m; j++) + { + for (unsigned int i = 0; i < n; i++) + deallog << y[j].getADValue(i) << (i -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -std::ofstream logfile("output"); - -template -void test (unsigned n_refinements) -{ - unsigned int rank = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - - if (rank == 0) - { - deallog << "Checking in " << dim << " space dimensions" - << std::endl - << "---------------------------------------" << std::endl; - } - - // create serial grid and refine once - parallel::distributed::Triangulation tria1(MPI_COMM_SELF); - GridGenerator::hyper_cube (tria1, -1, 1); - tria1.refine_global(1); - - // create distributed grid and refine twice - parallel::distributed::Triangulation tria2(MPI_COMM_WORLD); - GridGenerator::hyper_cube (tria2, -1, 1); - tria2.refine_global(n_refinements); - - // do some local refinement - Point p0; - p0 *= 0.; - for (int i = 0; i < n_refinements; ++i) - { - typename Triangulation::active_cell_iterator cell; - for (cell = tria2.begin_active(); cell != tria2.end(); ++cell) - { - if (cell->is_locally_owned() && (cell->center().distance(p0) < 0.71 / double(i + 1))) - cell->set_refine_flag(); - } - - tria2.prepare_coarsening_and_refinement (); - tria2.execute_coarsening_and_refinement (); - } - - FE_DGQ fe_dg(0); - DoFHandler dof_handler1(tria1); - dof_handler1.distribute_dofs(fe_dg); - - DoFHandler dof_handler2(tria2); - dof_handler2.distribute_dofs(fe_dg); - - InterGridMap > grid_1_to_2_map; - grid_1_to_2_map.make_mapping (dof_handler1, dof_handler2); - - typedef std::vector > TransferRep; - TransferRep transfer_representation; - DoFTools::compute_intergrid_transfer_representation (dof_handler1, 0, dof_handler2, - 0, grid_1_to_2_map, - transfer_representation); - - // For this test case, all weights are one and their sum - // should be equal to number of degrees of freedom - unsigned local_sum = 0.; - for (size_t i = 0; i < transfer_representation.size(); ++i) - { - TransferRep::value_type m = transfer_representation[i]; - for (TransferRep::value_type::const_iterator it = m.begin(); it != m.end(); ++it) - local_sum += it->second; - } - - unsigned global_sum = Utilities::MPI::sum(local_sum, MPI_COMM_WORLD); - - if (rank == 0) - { - deallog << "# dofs = " << dof_handler2.n_dofs() << std::endl; - deallog << "sum(weights) = " << global_sum << std::endl; - if (dof_handler2.n_dofs() == global_sum) - deallog << "OK" << std::endl; - else - deallog << "Failed" << std::endl; - } -} - -int main (int argc, char *argv[]) -{ - Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); - - unsigned int rank = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - - deallog.push(Utilities::int_to_string(rank)); - - unsigned n_refinements = 2; - - if (rank == 0) - { - std::ofstream logfile("output"); - deallog.attach(logfile); - - deallog.push("2d"); - test<2>(n_refinements); - deallog.pop(); - deallog.push("3d"); - test<3>(n_refinements); - deallog.pop(); - } - else - { - test<2>(n_refinements); - test<3>(n_refinements); - } -} +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Purpose: checks compute_intergrid_transfer_representation + * on distributed (p4est) triangulation + * + * Author: Alexander Grayver, 2015 + */ + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +std::ofstream logfile("output"); + +template +void test (unsigned n_refinements) +{ + unsigned int rank = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + if (rank == 0) + { + deallog << "Checking in " << dim << " space dimensions" + << std::endl + << "---------------------------------------" << std::endl; + } + + // create serial grid and refine once + parallel::distributed::Triangulation tria1(MPI_COMM_SELF); + GridGenerator::hyper_cube (tria1, -1, 1); + tria1.refine_global(1); + + // create distributed grid and refine twice + parallel::distributed::Triangulation tria2(MPI_COMM_WORLD); + GridGenerator::hyper_cube (tria2, -1, 1); + tria2.refine_global(n_refinements); + + // do some local refinement + Point p0; + p0 *= 0.; + for (int i = 0; i < n_refinements; ++i) + { + typename Triangulation::active_cell_iterator cell; + for (cell = tria2.begin_active(); cell != tria2.end(); ++cell) + { + if (cell->is_locally_owned() && (cell->center().distance(p0) < 0.71 / double(i + 1))) + cell->set_refine_flag(); + } + + tria2.prepare_coarsening_and_refinement (); + tria2.execute_coarsening_and_refinement (); + } + + FE_DGQ fe_dg(0); + DoFHandler dof_handler1(tria1); + dof_handler1.distribute_dofs(fe_dg); + + DoFHandler dof_handler2(tria2); + dof_handler2.distribute_dofs(fe_dg); + + InterGridMap > grid_1_to_2_map; + grid_1_to_2_map.make_mapping (dof_handler1, dof_handler2); + + typedef std::vector > TransferRep; + TransferRep transfer_representation; + DoFTools::compute_intergrid_transfer_representation (dof_handler1, 0, dof_handler2, + 0, grid_1_to_2_map, + transfer_representation); + + // For this test case, all weights are one and their sum + // should be equal to number of degrees of freedom + unsigned local_sum = 0.; + for (size_t i = 0; i < transfer_representation.size(); ++i) + { + TransferRep::value_type m = transfer_representation[i]; + for (TransferRep::value_type::const_iterator it = m.begin(); it != m.end(); ++it) + local_sum += it->second; + } + + unsigned global_sum = Utilities::MPI::sum(local_sum, MPI_COMM_WORLD); + + if (rank == 0) + { + deallog << "# dofs = " << dof_handler2.n_dofs() << std::endl; + deallog << "sum(weights) = " << global_sum << std::endl; + if (dof_handler2.n_dofs() == global_sum) + deallog << "OK" << std::endl; + else + deallog << "Failed" << std::endl; + } +} + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + unsigned int rank = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + deallog.push(Utilities::int_to_string(rank)); + + unsigned n_refinements = 2; + + if (rank == 0) + { + std::ofstream logfile("output"); + deallog.attach(logfile); + + deallog.push("2d"); + test<2>(n_refinements); + deallog.pop(); + deallog.push("3d"); + test<3>(n_refinements); + deallog.pop(); + } + else + { + test<2>(n_refinements); + test<3>(n_refinements); + } +} diff --git a/tests/fe/fe_project_2d.cc b/tests/fe/fe_project_2d.cc index 7d0b15c021..696b68a923 100644 --- a/tests/fe/fe_project_2d.cc +++ b/tests/fe/fe_project_2d.cc @@ -1,263 +1,263 @@ -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../tests.h" - -/* - * This program projects a function into FE spaces defined on meshes - * consisting of: rectangular cells, affine cells, non-affine cells - * - * The error, curl and divergence are then numerically calculated on a - * series of globally refined meshes and get output. - * - * Among FE spaces tested are: FE_ABF, FE_Nedelec, FE_RaviartThomas, - * FE_Q^dim (via FESystem) - * - * Alexander Grayver - */ - -using namespace dealii; - -static const Point<2> vertices_nonaffine[] = -{ - Point<2> (-1., -1.), - Point<2> (0., -1.), - Point<2> (1., -1.), - - Point<2> (-1., 0.), - Point<2> (0.3, 0.3), - Point<2> (1., 0.), - - Point<2> (-1., 1.), - Point<2> (0., 1.), - Point<2> (1., 1.), -}; - -static const Point<2> vertices_affine[] = -{ - Point<2> (-1.4, -1.), - Point<2> (-0.4, -1.), - Point<2> (0.6, -1.), - - Point<2> (-1.2, 0.), - Point<2> (-0.2, 0.), - Point<2> (0.8, 0.), - - Point<2> (-1., 1.), - Point<2> (0., 1.), - Point<2> (1., 1.), -}; - -static const Point<2> vertices_rectangular[] = -{ - Point<2> (-1., -1.), - Point<2> (0., -1.), - Point<2> (1., -1.), - - Point<2> (-1., 0.), - Point<2> (0., 0.), - Point<2> (1., 0.), - - Point<2> (-1., 1.), - Point<2> (0., 1.), - Point<2> (1., 1.), -}; - -static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]); -static const unsigned n_cells = 4; - -template -class VectorFunction : public Function -{ -public: - VectorFunction() : Function(dim) {} - virtual double value (const Point &p, const unsigned int component) const; - virtual void vector_value(const Point &p, Vector &values) const; -}; - -template <> -double VectorFunction<2>::value(const Point<2> &p, const unsigned int component) const -{ - Assert (component < 2, ExcIndexRange (component, 0, 1)); - - const double PI = numbers::PI; - double val = 0.0; - switch (component) - { - case 0: - val = cos(PI*p(0))*sin(PI*p(1)); - break; - case 1: - val = -sin(PI*p(0))*cos(PI*p(1)); - break; - } - return val; -} - -template -void VectorFunction::vector_value(const Point &p, Vector &values) const -{ - for (int i = 0; i < dim; ++i) - values(i) = value(p, i); -} - -void create_tria(Triangulation<2> &triangulation, const Point<2> *vertices_parallelograms) -{ - const std::vector > vertices (&vertices_parallelograms[0], - &vertices_parallelograms[n_vertices]); - - static const int cell_vertices[][GeometryInfo<2>::vertices_per_cell] = - { - {0, 1, 3, 4}, - {1, 2, 4, 5}, - {3, 4, 6, 7}, - {4, 5, 7, 8} - }; - - std::vector > cells (n_cells, CellData<2>()); - for (unsigned i = 0; i::vertices_per_cell; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - } - - triangulation.create_triangulation (vertices, cells, SubCellData()); - triangulation.refine_global(1); -} - -template -void test(const FiniteElement &fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) -{ - deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; - deallog << "DoFs\t\t||u-u_h||\tcurl(u_h)\tdiv(u_h)" << std::endl; - - Triangulation triangulation; - create_tria(triangulation, vertices_parallelograms); - - DoFHandler dof_handler(triangulation); - - VectorFunction fe_function; - const FEValuesExtractors::Vector vec (0); - const QGauss quadrature (fe.degree+1); - const unsigned int n_q_points = quadrature.size (); - MappingQ mapping(1); - //MappingQGeneric mapping(1); - std::vector div_v(n_q_points); - std::vector::curl_type> curl_v(n_q_points); - - for (unsigned cycle = 0; cycle < n_cycles; ++cycle) - { - dof_handler.distribute_dofs(fe); - - ConstraintMatrix constraints; - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - constraints.close(); - - Vector v(dof_handler.n_dofs()); - VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); - - Vector diff(triangulation.n_active_cells()); - VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, - QGauss(fe.degree + 2), VectorTools::L2_norm); - - typename FEValuesViews::Vector::curl_type total_curl; - double total_div = 0; - total_curl *= 0.; - - FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients); - unsigned int cell_index = 0; - - for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); - cell != dof_handler.end (); ++cell, ++cell_index) - { - fe_values.reinit (cell); - const std::vector &JxW_values = fe_values.get_JxW_values (); - fe_values[vec].get_function_divergences (v, div_v); - fe_values[vec].get_function_curls (v, curl_v); - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - { - total_div += JxW_values[q_point] * div_v[q_point]; - total_curl += JxW_values[q_point] * curl_v[q_point]; - } - } - - deallog << dof_handler.n_dofs() << "\t\t" - << diff.l2_norm() << "\t" - << total_curl << "\t" - << total_div << std::endl; - - if (global) - triangulation.refine_global(); - else - { - GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); - triangulation.prepare_coarsening_and_refinement(); - triangulation.execute_coarsening_and_refinement(); - } - } -} - -int main () -{ - std::ofstream logfile ("output"); - deallog << std::setprecision(7); - deallog << std::fixed; - deallog.attach(logfile); - - const static unsigned dim = 2; - unsigned order = 1; - unsigned n_cycles = 4; - - deallog << "2d\nRectangular grid:\n"; - - const Point *vertices = &vertices_rectangular[0]; - test(FE_Nedelec (order), n_cycles, true, vertices); - test(FE_RaviartThomas (order), n_cycles, true, vertices); - test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); - test(FE_ABF (order), n_cycles, true, vertices); - - deallog << "\nAffine grid:\n"; - - vertices = &vertices_affine[0]; - test(FE_Nedelec (order), n_cycles, true, vertices); - test(FE_RaviartThomas (order), n_cycles, true, vertices); - test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); - test(FE_ABF (order), n_cycles, true, vertices); - - deallog << "\nNon-affine grid:\n"; - - vertices = &vertices_nonaffine[0]; - test(FE_Nedelec (order), n_cycles, true, vertices); - test(FE_RaviartThomas (order), n_cycles, true, vertices); - test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); - test(FE_ABF (order), n_cycles, true, vertices); - - deallog << std::endl; -} +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +/* + * This program projects a function into FE spaces defined on meshes + * consisting of: rectangular cells, affine cells, non-affine cells + * + * The error, curl and divergence are then numerically calculated on a + * series of globally refined meshes and get output. + * + * Among FE spaces tested are: FE_ABF, FE_Nedelec, FE_RaviartThomas, + * FE_Q^dim (via FESystem) + * + * Alexander Grayver + */ + +using namespace dealii; + +static const Point<2> vertices_nonaffine[] = +{ + Point<2> (-1., -1.), + Point<2> (0., -1.), + Point<2> (1., -1.), + + Point<2> (-1., 0.), + Point<2> (0.3, 0.3), + Point<2> (1., 0.), + + Point<2> (-1., 1.), + Point<2> (0., 1.), + Point<2> (1., 1.), +}; + +static const Point<2> vertices_affine[] = +{ + Point<2> (-1.4, -1.), + Point<2> (-0.4, -1.), + Point<2> (0.6, -1.), + + Point<2> (-1.2, 0.), + Point<2> (-0.2, 0.), + Point<2> (0.8, 0.), + + Point<2> (-1., 1.), + Point<2> (0., 1.), + Point<2> (1., 1.), +}; + +static const Point<2> vertices_rectangular[] = +{ + Point<2> (-1., -1.), + Point<2> (0., -1.), + Point<2> (1., -1.), + + Point<2> (-1., 0.), + Point<2> (0., 0.), + Point<2> (1., 0.), + + Point<2> (-1., 1.), + Point<2> (0., 1.), + Point<2> (1., 1.), +}; + +static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]); +static const unsigned n_cells = 4; + +template +class VectorFunction : public Function +{ +public: + VectorFunction() : Function(dim) {} + virtual double value (const Point &p, const unsigned int component) const; + virtual void vector_value(const Point &p, Vector &values) const; +}; + +template <> +double VectorFunction<2>::value(const Point<2> &p, const unsigned int component) const +{ + Assert (component < 2, ExcIndexRange (component, 0, 1)); + + const double PI = numbers::PI; + double val = 0.0; + switch (component) + { + case 0: + val = cos(PI*p(0))*sin(PI*p(1)); + break; + case 1: + val = -sin(PI*p(0))*cos(PI*p(1)); + break; + } + return val; +} + +template +void VectorFunction::vector_value(const Point &p, Vector &values) const +{ + for (int i = 0; i < dim; ++i) + values(i) = value(p, i); +} + +void create_tria(Triangulation<2> &triangulation, const Point<2> *vertices_parallelograms) +{ + const std::vector > vertices (&vertices_parallelograms[0], + &vertices_parallelograms[n_vertices]); + + static const int cell_vertices[][GeometryInfo<2>::vertices_per_cell] = + { + {0, 1, 3, 4}, + {1, 2, 4, 5}, + {3, 4, 6, 7}, + {4, 5, 7, 8} + }; + + std::vector > cells (n_cells, CellData<2>()); + for (unsigned i = 0; i::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation (vertices, cells, SubCellData()); + triangulation.refine_global(1); +} + +template +void test(const FiniteElement &fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) +{ + deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; + deallog << "DoFs\t\t||u-u_h||\tcurl(u_h)\tdiv(u_h)" << std::endl; + + Triangulation triangulation; + create_tria(triangulation, vertices_parallelograms); + + DoFHandler dof_handler(triangulation); + + VectorFunction fe_function; + const FEValuesExtractors::Vector vec (0); + const QGauss quadrature (fe.degree+1); + const unsigned int n_q_points = quadrature.size (); + MappingQ mapping(1); + //MappingQGeneric mapping(1); + std::vector div_v(n_q_points); + std::vector::curl_type> curl_v(n_q_points); + + for (unsigned cycle = 0; cycle < n_cycles; ++cycle) + { + dof_handler.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + Vector v(dof_handler.n_dofs()); + VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); + + Vector diff(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, + QGauss(fe.degree + 2), VectorTools::L2_norm); + + typename FEValuesViews::Vector::curl_type total_curl; + double total_div = 0; + total_curl *= 0.; + + FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients); + unsigned int cell_index = 0; + + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + cell != dof_handler.end (); ++cell, ++cell_index) + { + fe_values.reinit (cell); + const std::vector &JxW_values = fe_values.get_JxW_values (); + fe_values[vec].get_function_divergences (v, div_v); + fe_values[vec].get_function_curls (v, curl_v); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + total_div += JxW_values[q_point] * div_v[q_point]; + total_curl += JxW_values[q_point] * curl_v[q_point]; + } + } + + deallog << dof_handler.n_dofs() << "\t\t" + << diff.l2_norm() << "\t" + << total_curl << "\t" + << total_div << std::endl; + + if (global) + triangulation.refine_global(); + else + { + GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement(); + } + } +} + +int main () +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(7); + deallog << std::fixed; + deallog.attach(logfile); + + const static unsigned dim = 2; + unsigned order = 1; + unsigned n_cycles = 4; + + deallog << "2d\nRectangular grid:\n"; + + const Point *vertices = &vertices_rectangular[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); + test(FE_ABF (order), n_cycles, true, vertices); + + deallog << "\nAffine grid:\n"; + + vertices = &vertices_affine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); + test(FE_ABF (order), n_cycles, true, vertices); + + deallog << "\nNon-affine grid:\n"; + + vertices = &vertices_nonaffine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); + test(FE_ABF (order), n_cycles, true, vertices); + + deallog << std::endl; +} diff --git a/tests/fe/fe_project_3d.cc b/tests/fe/fe_project_3d.cc index 590b25d486..864b212be2 100644 --- a/tests/fe/fe_project_3d.cc +++ b/tests/fe/fe_project_3d.cc @@ -1,429 +1,429 @@ -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../tests.h" - -/* - * This program projects a function into FE spaces defined on meshes - * consisting of: rectangular cells, affine cells, non-affine cells - * - * The error, curl and divergence are then numerically calculated - * on a series of globally refined meshes and output. - * - * Among FE spaces tested are FE_Nedelec and FE_RaviartThomas - * - * Alexander Grayver, Maien Hamed - */ - -using namespace dealii; - -static const Point<3> vertices_affine[] = -{ - Point<3> (-1., -1., -1.), - Point<3> (0., -1., -1.), - Point<3> (1., -1., -1.), - - Point<3> (-1.2, -1., 0.), - Point<3> (-0.2, -1., 0.), - Point<3> (0.8, -1., 0.), - - Point<3> (-1.4, -1., 1), - Point<3> (-0.4, -1., 1), - Point<3> (0.6, -1., 1), - - Point<3> (-1., 0., -1.), - Point<3> (0., 0., -1.), - Point<3> (1., 0., -1.), - - Point<3> (-1.2, 0., 0.), - Point<3> (-0.2, 0., 0.), - Point<3> (0.8, 0., 0.), - - Point<3> (-1.4, 0., 1), - Point<3> (-0.4, 0., 1), - Point<3> (0.6, 0., 1), - - Point<3> (-1., 1., -1.), - Point<3> (0., 1., -1.), - Point<3> (1., 1., -1.), - - Point<3> (-1.2, 1., 0.), - Point<3> (-0.2, 1., 0.), - Point<3> (0.8, 1., 0.), - - Point<3> (-1.4, 1., 1), - Point<3> (-0.4, 1., 1), - Point<3> (0.6, 1., 1) -}; - -static const Point<3> vertices_nonaffine[] = -{ - Point<3> (-1., -1., -1.), - Point<3> (0., -1., -1.), - Point<3> (1., -1., -1.), - - Point<3> (-1., -1., 0.), - Point<3> (0., -1., 0.), - Point<3> (1., -1., 0.), - - Point<3> (-1., -1., 1), - Point<3> (0., -1., 1), - Point<3> (1., -1., 1), - - Point<3> (-1., 0., -1.), - Point<3> (0., 0., -1.), - Point<3> (1., 0., -1.), - - Point<3> (-1., 0., 0.), - Point<3> (0.2, 0.3, 0.1), - Point<3> (1., 0., 0.), - - Point<3> (-1., 0., 1), - Point<3> (0., 0., 1), - Point<3> (1., 0., 1), - - Point<3> (-1., 1., -1.), - Point<3> (0., 1., -1.), - Point<3> (1., 1., -1.), - - Point<3> (-1., 1., 0.), - Point<3> (0., 1., 0.), - Point<3> (1., 1., 0.), - - Point<3> (-1., 1., 1.), - Point<3> (0., 1., 1.), - Point<3> (1., 1., 1.) -}; - -static const Point<3> vertices_rectangular[] = -{ - Point<3> (-1., -1., -1.), - Point<3> (0., -1., -1.), - Point<3> (1., -1., -1.), - - Point<3> (-1., -1., 0.), - Point<3> (0., -1., 0.), - Point<3> (1., -1., 0.), - - Point<3> (-1., -1., 1), - Point<3> (0., -1., 1), - Point<3> (1., -1., 1), - - Point<3> (-1., 0., -1.), - Point<3> (0., 0., -1.), - Point<3> (1., 0., -1.), - - Point<3> (-1., 0., 0.), - Point<3> (0., 0., 0.), - Point<3> (1., 0., 0.), - - Point<3> (-1., 0., 1), - Point<3> (0., 0., 1), - Point<3> (1., 0., 1), - - Point<3> (-1., 1., -1.), - Point<3> (0., 1., -1.), - Point<3> (1., 1., -1.), - - Point<3> (-1., 1., 0.), - Point<3> (0., 1., 0.), - Point<3> (1., 1., 0.), - - Point<3> (-1., 1., 1.), - Point<3> (0., 1., 1.), - Point<3> (1., 1., 1.) -}; - -static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]); -static const unsigned n_cells = 8; - -template -class VectorFunction : public Function -{ -public: - VectorFunction() : Function(dim) {} - virtual double value (const Point &p, const unsigned int component) const; - virtual void vector_value(const Point &p, Vector &values) const; - virtual Tensor< 1, dim > gradient (const Point< dim > &p, const unsigned int component=0) const; -}; - -template <> -double VectorFunction<3>::value(const Point<3> &p, const unsigned int component) const -{ - Assert (component < 3, ExcIndexRange (component, 0, 2)); - - const double PI = numbers::PI; - double val = 0.0; - switch (component) - { - case 0: - val = -sin(PI*p(0))*cos(PI*p(1))*cos(PI*p(2)); - break; - case 1: - val = -cos(PI*p(0))*sin(PI*p(1))*cos(PI*p(2)); - break; - case 2: - val = 2*cos(PI*p(0))*cos(PI*p(1))*sin(PI*p(2)); - break; - } - return val; -} - -template -void VectorFunction::vector_value(const Point &p, Vector &values) const -{ - for (int i = 0; i < dim; ++i) - values(i) = value(p, i); -} - -template <> -Tensor<1, 3> VectorFunction<3>::gradient(const Point<3> &p, const unsigned int component) const -{ - const double PI = numbers::PI; - Tensor<1, 3> val; - double x = p(0), y = p(1), z = p(2); - - switch (component) - { - case 0: - val[0] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); - val[1] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); - val[2] = -2*PI*cos(PI*y)*sin(PI*x)*sin(PI*z); - break; - case 1: - val[0] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); - val[1] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); - val[2] = -2*PI*cos(PI*x)*sin(PI*y)*sin(PI*z); - break; - case 2: - val[0] = PI*cos(PI*y)*sin(PI*x)*sin(PI*z); - val[1] = PI*cos(PI*x)*sin(PI*y)*sin(PI*z); - val[2] = 2*PI*cos(PI*x)*cos(PI*y)*cos(PI*z); - break; - } - return val; -} - -void create_tria(Triangulation<3> &triangulation, const Point<3> *vertices_parallelograms) -{ - const std::vector > vertices (&vertices_parallelograms[0], - &vertices_parallelograms[n_vertices]); - - // create grid with all possible combintations of face_flip, face_orientation and face_rotation flags - static const int cell_vertices[][GeometryInfo<3>::vertices_per_cell] = - { - {0, 1, 9, 10, 3, 4, 12, 13}, // cell 1 standard - {1, 2, 10, 11, 4, 5, 13, 14}, // cell 2 standard - //{10, 11, 13, 14, 1, 2, 4, 5}, // cell 2 rotated by 270 deg - {9, 10, 18, 19, 12, 13, 21, 22}, // cell 3 standard - {10, 11, 19, 20, 13, 14, 22, 23}, // cell 4 standard - //{13, 14, 10, 11, 22, 23, 19, 20}, // cell 4 rotated by 90 deg - {3, 4, 12, 13, 6, 7, 15, 16}, // cell 5 standard - {4, 5, 13, 14, 7, 8, 16, 17}, // cell 6 standard - {12, 13, 21, 22, 15, 16, 24, 25}, // cell 7 standard - //{24, 25, 15, 16, 21, 22, 12, 13}, // cell 7 rotated by 180 deg - {13, 14, 22, 23, 16, 17, 25, 26} // cell 8 standard - }; - - std::vector > cells (n_cells, CellData<3>()); - for (unsigned i = 0; i::vertices_per_cell; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - } - - triangulation.create_triangulation (vertices, cells, SubCellData()); -} - -template -void test(const FiniteElement &fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) -{ - deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; - deallog << "DoFs\t\t||u-u_h||_1\tcurl(u_h)\ttangentials\tcurl(curl(u_h))\tcurl_curl_traces\tdiv(u_h)\tboundary_flux" << std::endl; - - Triangulation triangulation; - create_tria(triangulation, vertices_parallelograms); - - DoFHandler dof_handler(triangulation); - - VectorFunction fe_function; - const FEValuesExtractors::Vector vec (0); - const QGauss quadrature (fe.degree+1); - const QGauss face_quadrature (fe.degree+1); - const unsigned int n_q_points = quadrature.size (); - const unsigned int n_face_q_points = face_quadrature.size (); - //MappingQ mapping(2); - MappingQGeneric mapping(1); - std::vector div_v(n_q_points); - std::vector::curl_type> curl_v(n_q_points); - std::vector > hessians(n_q_points); - - std::vector > face_values (n_face_q_points); - std::vector::curl_type> face_curls (n_face_q_points); - - for (unsigned cycle = 0; cycle < n_cycles; ++cycle) - { - dof_handler.distribute_dofs(fe); - - ConstraintMatrix constraints; - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - constraints.close(); - - Vector v(dof_handler.n_dofs()); - VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); - - Vector diff(triangulation.n_active_cells()); - VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, - QGauss(fe.degree + 2), VectorTools::L1_norm); - - typename FEValuesViews::Vector::curl_type total_curl, boundary_tangentials; - Tensor<1, dim> total_curl_curl, boundary_curl_curl_traces; - double total_div = 0; - double boundary_flux = 0; - total_curl *= 0.; - boundary_tangentials *= 0.; - - FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_hessians); - FEFaceValues fe_face_values(mapping, fe, face_quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_normal_vectors ); - unsigned int cell_index = 0; - - for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); - cell != dof_handler.end (); ++cell, ++cell_index) - { - fe_values.reinit (cell); - const std::vector &JxW_values = fe_values.get_JxW_values (); - fe_values[vec].get_function_divergences (v, div_v); - fe_values[vec].get_function_curls (v, curl_v); - fe_values[vec].get_function_hessians (v, hessians); - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - { - total_div += JxW_values[q_point] * div_v[q_point]; - total_curl += JxW_values[q_point] * curl_v[q_point]; - if (dim == 3) - { - total_curl_curl[0] += JxW_values[q_point] * (hessians[q_point][1][0][1] + hessians[q_point][2][0][2] - hessians[q_point][0][1][1] - hessians[q_point][0][2][2]); - total_curl_curl[1] += JxW_values[q_point] * (hessians[q_point][2][1][2] + hessians[q_point][0][0][1] - hessians[q_point][1][2][2] - hessians[q_point][1][0][0]); - total_curl_curl[2] += JxW_values[q_point] * (hessians[q_point][0][0][2] + hessians[q_point][1][1][2] - hessians[q_point][2][0][0] - hessians[q_point][2][1][1]); - } - } - - for (unsigned int face=0; face::faces_per_cell; ++face) - { - fe_face_values.reinit(cell,face); - const std::vector &face_JxW_values = fe_face_values.get_JxW_values (); - fe_face_values[vec].get_function_values (v, face_values); - if (dim==3) fe_face_values[vec].get_function_curls (v, face_curls); - for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) - { - const Tensor<1,dim> &normal = fe_face_values.normal_vector(q_point); - - // boundary flux - if (cell->at_boundary(face)) - boundary_flux += face_JxW_values[q_point] * (face_values[q_point] * normal); - else - total_div -= face_JxW_values[q_point] * (face_values[q_point] * normal); - - // boundary tangentials (curl traces) - typename FEValuesViews::Vector::curl_type n_x_v; - if (dim==2) - n_x_v[0] = (-normal[1]*face_values[q_point][0] + normal[0]*face_values[q_point][1]); - else if (dim==3) - cross_product(*reinterpret_cast*>(&n_x_v), normal, face_values[q_point]); - - if (cell->at_boundary(face)) - boundary_tangentials += face_JxW_values[q_point] * n_x_v; - else - total_curl -= face_JxW_values[q_point] * n_x_v; - - // boundary curl curl traces - if (dim==3) - { - Tensor<1,dim> n_x_curl_u; - cross_product(n_x_curl_u, normal, *reinterpret_cast*>(&face_curls[q_point])); - if (cell->at_boundary(face)) - boundary_curl_curl_traces += face_JxW_values[q_point] * n_x_curl_u; - else - total_curl_curl -= face_JxW_values[q_point] * n_x_curl_u; - } - } - } - } - - deallog << dof_handler.n_dofs() << "\t\t" - << diff.l1_norm() << "\t" - << total_curl.norm() << "\t" - << boundary_tangentials.norm() << "\t" - << total_curl_curl.norm() << "\t" - << boundary_curl_curl_traces.norm() << "\t" - << total_div << "\t" - << boundary_flux << std::endl; - - if (global) - triangulation.refine_global(); - else - { - GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); - triangulation.prepare_coarsening_and_refinement(); - triangulation.execute_coarsening_and_refinement(); - } - } -} - -int main () -{ - std::ofstream logfile ("output"); - deallog << std::setprecision(7); - deallog << std::fixed; - deallog.attach(logfile); - - const static unsigned dim = 3; - unsigned order = 1; - unsigned n_cycles = 2; - - deallog << "3d\nRectangular grid:\n"; - - const Point *vertices = &vertices_rectangular[0]; - test(FE_Nedelec (order), n_cycles, true, vertices); - test(FE_RaviartThomas (order), n_cycles, true, vertices); - - deallog << "\nAffine grid:\n"; - - vertices = &vertices_affine[0]; - test(FE_Nedelec (order), n_cycles, true, vertices); - test(FE_RaviartThomas (order), n_cycles, true, vertices); - - deallog << "\nNon-affine grid:\n"; - - vertices = &vertices_nonaffine[0]; - test(FE_Nedelec (order), n_cycles, true, vertices); - test(FE_RaviartThomas (order), n_cycles, true, vertices); - - deallog << std::endl; -} +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +/* + * This program projects a function into FE spaces defined on meshes + * consisting of: rectangular cells, affine cells, non-affine cells + * + * The error, curl and divergence are then numerically calculated + * on a series of globally refined meshes and output. + * + * Among FE spaces tested are FE_Nedelec and FE_RaviartThomas + * + * Alexander Grayver, Maien Hamed + */ + +using namespace dealii; + +static const Point<3> vertices_affine[] = +{ + Point<3> (-1., -1., -1.), + Point<3> (0., -1., -1.), + Point<3> (1., -1., -1.), + + Point<3> (-1.2, -1., 0.), + Point<3> (-0.2, -1., 0.), + Point<3> (0.8, -1., 0.), + + Point<3> (-1.4, -1., 1), + Point<3> (-0.4, -1., 1), + Point<3> (0.6, -1., 1), + + Point<3> (-1., 0., -1.), + Point<3> (0., 0., -1.), + Point<3> (1., 0., -1.), + + Point<3> (-1.2, 0., 0.), + Point<3> (-0.2, 0., 0.), + Point<3> (0.8, 0., 0.), + + Point<3> (-1.4, 0., 1), + Point<3> (-0.4, 0., 1), + Point<3> (0.6, 0., 1), + + Point<3> (-1., 1., -1.), + Point<3> (0., 1., -1.), + Point<3> (1., 1., -1.), + + Point<3> (-1.2, 1., 0.), + Point<3> (-0.2, 1., 0.), + Point<3> (0.8, 1., 0.), + + Point<3> (-1.4, 1., 1), + Point<3> (-0.4, 1., 1), + Point<3> (0.6, 1., 1) +}; + +static const Point<3> vertices_nonaffine[] = +{ + Point<3> (-1., -1., -1.), + Point<3> (0., -1., -1.), + Point<3> (1., -1., -1.), + + Point<3> (-1., -1., 0.), + Point<3> (0., -1., 0.), + Point<3> (1., -1., 0.), + + Point<3> (-1., -1., 1), + Point<3> (0., -1., 1), + Point<3> (1., -1., 1), + + Point<3> (-1., 0., -1.), + Point<3> (0., 0., -1.), + Point<3> (1., 0., -1.), + + Point<3> (-1., 0., 0.), + Point<3> (0.2, 0.3, 0.1), + Point<3> (1., 0., 0.), + + Point<3> (-1., 0., 1), + Point<3> (0., 0., 1), + Point<3> (1., 0., 1), + + Point<3> (-1., 1., -1.), + Point<3> (0., 1., -1.), + Point<3> (1., 1., -1.), + + Point<3> (-1., 1., 0.), + Point<3> (0., 1., 0.), + Point<3> (1., 1., 0.), + + Point<3> (-1., 1., 1.), + Point<3> (0., 1., 1.), + Point<3> (1., 1., 1.) +}; + +static const Point<3> vertices_rectangular[] = +{ + Point<3> (-1., -1., -1.), + Point<3> (0., -1., -1.), + Point<3> (1., -1., -1.), + + Point<3> (-1., -1., 0.), + Point<3> (0., -1., 0.), + Point<3> (1., -1., 0.), + + Point<3> (-1., -1., 1), + Point<3> (0., -1., 1), + Point<3> (1., -1., 1), + + Point<3> (-1., 0., -1.), + Point<3> (0., 0., -1.), + Point<3> (1., 0., -1.), + + Point<3> (-1., 0., 0.), + Point<3> (0., 0., 0.), + Point<3> (1., 0., 0.), + + Point<3> (-1., 0., 1), + Point<3> (0., 0., 1), + Point<3> (1., 0., 1), + + Point<3> (-1., 1., -1.), + Point<3> (0., 1., -1.), + Point<3> (1., 1., -1.), + + Point<3> (-1., 1., 0.), + Point<3> (0., 1., 0.), + Point<3> (1., 1., 0.), + + Point<3> (-1., 1., 1.), + Point<3> (0., 1., 1.), + Point<3> (1., 1., 1.) +}; + +static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]); +static const unsigned n_cells = 8; + +template +class VectorFunction : public Function +{ +public: + VectorFunction() : Function(dim) {} + virtual double value (const Point &p, const unsigned int component) const; + virtual void vector_value(const Point &p, Vector &values) const; + virtual Tensor< 1, dim > gradient (const Point< dim > &p, const unsigned int component=0) const; +}; + +template <> +double VectorFunction<3>::value(const Point<3> &p, const unsigned int component) const +{ + Assert (component < 3, ExcIndexRange (component, 0, 2)); + + const double PI = numbers::PI; + double val = 0.0; + switch (component) + { + case 0: + val = -sin(PI*p(0))*cos(PI*p(1))*cos(PI*p(2)); + break; + case 1: + val = -cos(PI*p(0))*sin(PI*p(1))*cos(PI*p(2)); + break; + case 2: + val = 2*cos(PI*p(0))*cos(PI*p(1))*sin(PI*p(2)); + break; + } + return val; +} + +template +void VectorFunction::vector_value(const Point &p, Vector &values) const +{ + for (int i = 0; i < dim; ++i) + values(i) = value(p, i); +} + +template <> +Tensor<1, 3> VectorFunction<3>::gradient(const Point<3> &p, const unsigned int component) const +{ + const double PI = numbers::PI; + Tensor<1, 3> val; + double x = p(0), y = p(1), z = p(2); + + switch (component) + { + case 0: + val[0] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + val[1] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); + val[2] = -2*PI*cos(PI*y)*sin(PI*x)*sin(PI*z); + break; + case 1: + val[0] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); + val[1] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + val[2] = -2*PI*cos(PI*x)*sin(PI*y)*sin(PI*z); + break; + case 2: + val[0] = PI*cos(PI*y)*sin(PI*x)*sin(PI*z); + val[1] = PI*cos(PI*x)*sin(PI*y)*sin(PI*z); + val[2] = 2*PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + break; + } + return val; +} + +void create_tria(Triangulation<3> &triangulation, const Point<3> *vertices_parallelograms) +{ + const std::vector > vertices (&vertices_parallelograms[0], + &vertices_parallelograms[n_vertices]); + + // create grid with all possible combintations of face_flip, face_orientation and face_rotation flags + static const int cell_vertices[][GeometryInfo<3>::vertices_per_cell] = + { + {0, 1, 9, 10, 3, 4, 12, 13}, // cell 1 standard + {1, 2, 10, 11, 4, 5, 13, 14}, // cell 2 standard + //{10, 11, 13, 14, 1, 2, 4, 5}, // cell 2 rotated by 270 deg + {9, 10, 18, 19, 12, 13, 21, 22}, // cell 3 standard + {10, 11, 19, 20, 13, 14, 22, 23}, // cell 4 standard + //{13, 14, 10, 11, 22, 23, 19, 20}, // cell 4 rotated by 90 deg + {3, 4, 12, 13, 6, 7, 15, 16}, // cell 5 standard + {4, 5, 13, 14, 7, 8, 16, 17}, // cell 6 standard + {12, 13, 21, 22, 15, 16, 24, 25}, // cell 7 standard + //{24, 25, 15, 16, 21, 22, 12, 13}, // cell 7 rotated by 180 deg + {13, 14, 22, 23, 16, 17, 25, 26} // cell 8 standard + }; + + std::vector > cells (n_cells, CellData<3>()); + for (unsigned i = 0; i::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation (vertices, cells, SubCellData()); +} + +template +void test(const FiniteElement &fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) +{ + deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; + deallog << "DoFs\t\t||u-u_h||_1\tcurl(u_h)\ttangentials\tcurl(curl(u_h))\tcurl_curl_traces\tdiv(u_h)\tboundary_flux" << std::endl; + + Triangulation triangulation; + create_tria(triangulation, vertices_parallelograms); + + DoFHandler dof_handler(triangulation); + + VectorFunction fe_function; + const FEValuesExtractors::Vector vec (0); + const QGauss quadrature (fe.degree+1); + const QGauss face_quadrature (fe.degree+1); + const unsigned int n_q_points = quadrature.size (); + const unsigned int n_face_q_points = face_quadrature.size (); + //MappingQ mapping(2); + MappingQGeneric mapping(1); + std::vector div_v(n_q_points); + std::vector::curl_type> curl_v(n_q_points); + std::vector > hessians(n_q_points); + + std::vector > face_values (n_face_q_points); + std::vector::curl_type> face_curls (n_face_q_points); + + for (unsigned cycle = 0; cycle < n_cycles; ++cycle) + { + dof_handler.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + Vector v(dof_handler.n_dofs()); + VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); + + Vector diff(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, + QGauss(fe.degree + 2), VectorTools::L1_norm); + + typename FEValuesViews::Vector::curl_type total_curl, boundary_tangentials; + Tensor<1, dim> total_curl_curl, boundary_curl_curl_traces; + double total_div = 0; + double boundary_flux = 0; + total_curl *= 0.; + boundary_tangentials *= 0.; + + FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_hessians); + FEFaceValues fe_face_values(mapping, fe, face_quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_normal_vectors ); + unsigned int cell_index = 0; + + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + cell != dof_handler.end (); ++cell, ++cell_index) + { + fe_values.reinit (cell); + const std::vector &JxW_values = fe_values.get_JxW_values (); + fe_values[vec].get_function_divergences (v, div_v); + fe_values[vec].get_function_curls (v, curl_v); + fe_values[vec].get_function_hessians (v, hessians); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + total_div += JxW_values[q_point] * div_v[q_point]; + total_curl += JxW_values[q_point] * curl_v[q_point]; + if (dim == 3) + { + total_curl_curl[0] += JxW_values[q_point] * (hessians[q_point][1][0][1] + hessians[q_point][2][0][2] - hessians[q_point][0][1][1] - hessians[q_point][0][2][2]); + total_curl_curl[1] += JxW_values[q_point] * (hessians[q_point][2][1][2] + hessians[q_point][0][0][1] - hessians[q_point][1][2][2] - hessians[q_point][1][0][0]); + total_curl_curl[2] += JxW_values[q_point] * (hessians[q_point][0][0][2] + hessians[q_point][1][1][2] - hessians[q_point][2][0][0] - hessians[q_point][2][1][1]); + } + } + + for (unsigned int face=0; face::faces_per_cell; ++face) + { + fe_face_values.reinit(cell,face); + const std::vector &face_JxW_values = fe_face_values.get_JxW_values (); + fe_face_values[vec].get_function_values (v, face_values); + if (dim==3) fe_face_values[vec].get_function_curls (v, face_curls); + for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) + { + const Tensor<1,dim> &normal = fe_face_values.normal_vector(q_point); + + // boundary flux + if (cell->at_boundary(face)) + boundary_flux += face_JxW_values[q_point] * (face_values[q_point] * normal); + else + total_div -= face_JxW_values[q_point] * (face_values[q_point] * normal); + + // boundary tangentials (curl traces) + typename FEValuesViews::Vector::curl_type n_x_v; + if (dim==2) + n_x_v[0] = (-normal[1]*face_values[q_point][0] + normal[0]*face_values[q_point][1]); + else if (dim==3) + cross_product(*reinterpret_cast*>(&n_x_v), normal, face_values[q_point]); + + if (cell->at_boundary(face)) + boundary_tangentials += face_JxW_values[q_point] * n_x_v; + else + total_curl -= face_JxW_values[q_point] * n_x_v; + + // boundary curl curl traces + if (dim==3) + { + Tensor<1,dim> n_x_curl_u; + cross_product(n_x_curl_u, normal, *reinterpret_cast*>(&face_curls[q_point])); + if (cell->at_boundary(face)) + boundary_curl_curl_traces += face_JxW_values[q_point] * n_x_curl_u; + else + total_curl_curl -= face_JxW_values[q_point] * n_x_curl_u; + } + } + } + } + + deallog << dof_handler.n_dofs() << "\t\t" + << diff.l1_norm() << "\t" + << total_curl.norm() << "\t" + << boundary_tangentials.norm() << "\t" + << total_curl_curl.norm() << "\t" + << boundary_curl_curl_traces.norm() << "\t" + << total_div << "\t" + << boundary_flux << std::endl; + + if (global) + triangulation.refine_global(); + else + { + GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement(); + } + } +} + +int main () +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(7); + deallog << std::fixed; + deallog.attach(logfile); + + const static unsigned dim = 3; + unsigned order = 1; + unsigned n_cycles = 2; + + deallog << "3d\nRectangular grid:\n"; + + const Point *vertices = &vertices_rectangular[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + + deallog << "\nAffine grid:\n"; + + vertices = &vertices_affine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + + deallog << "\nNon-affine grid:\n"; + + vertices = &vertices_nonaffine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + + deallog << std::endl; +} diff --git a/tests/numerics/no_flux_13.cc b/tests/numerics/no_flux_13.cc index 0f472b4dba..6e32dcdcb7 100644 --- a/tests/numerics/no_flux_13.cc +++ b/tests/numerics/no_flux_13.cc @@ -1,71 +1,71 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2012 - 2016 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -// bug report from mailing list from 11/15/2013 (simplified). no_normal_flux throws -// an ExcInternalError when handing it an FE with less than dim components. This -// is now fixed (throws ExcMessage). - -#include "../tests.h" - -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace dealii; - -template -void test() -{ - Triangulation triangulation; - GridGenerator::hyper_cube (triangulation,-1.0,1.0); - triangulation.begin_active()->face(1)->set_all_boundary_ids(1); - - FE_Q fe(1); - DoFHandler dof_handler(triangulation); - dof_handler.distribute_dofs(fe); - - ConstraintMatrix constraints; - std::set no_normal_flux_boundaries; - no_normal_flux_boundaries.insert (1); - deal_II_exceptions::disable_abort_on_exception(); - try - { - - VectorTools::compute_no_normal_flux_constraints (dof_handler, - 0, - no_normal_flux_boundaries, - constraints); - } - catch (ExceptionBase &e) - { - deallog << e.get_exc_name() << std::endl; - } - - constraints.close(); - constraints.print(deallog.get_file_stream ()); -} - - -int main () -{ - initlog(); - - test<3> (); -} +// --------------------------------------------------------------------- +// +// Copyright (C) 2012 - 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// bug report from mailing list from 11/15/2013 (simplified). no_normal_flux throws +// an ExcInternalError when handing it an FE with less than dim components. This +// is now fixed (throws ExcMessage). + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +template +void test() +{ + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation,-1.0,1.0); + triangulation.begin_active()->face(1)->set_all_boundary_ids(1); + + FE_Q fe(1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + ConstraintMatrix constraints; + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (1); + deal_II_exceptions::disable_abort_on_exception(); + try + { + + VectorTools::compute_no_normal_flux_constraints (dof_handler, + 0, + no_normal_flux_boundaries, + constraints); + } + catch (ExceptionBase &e) + { + deallog << e.get_exc_name() << std::endl; + } + + constraints.close(); + constraints.print(deallog.get_file_stream ()); +} + + +int main () +{ + initlog(); + + test<3> (); +} diff --git a/tests/numerics/project_bv_curl_conf_02.cc b/tests/numerics/project_bv_curl_conf_02.cc index 118514f036..46552e0737 100644 --- a/tests/numerics/project_bv_curl_conf_02.cc +++ b/tests/numerics/project_bv_curl_conf_02.cc @@ -1,130 +1,130 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2014 - 2015 by the Alexander Grayver & deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -// The test checks that project_boundary_values_curl_conforming -// works correctly for high-order FE_Nedelec elements used via -// FESystem. This requires the produced constraints to be the same -// for FE_Nedelec and FESystem(FE_Nedelec, 1). - -#include "../tests.h" -#include -#include -#include -#include -#include -#include -#include -#include - -std::ofstream logfile ("output"); - -template -class BoundaryFunction: public Function -{ -public: - BoundaryFunction (); - virtual void vector_value (const Point &p, Vector &values) const; -}; - -template -BoundaryFunction::BoundaryFunction (): Function (dim) -{ -} - -template -void BoundaryFunction::vector_value (const Point &, - Vector &values) const -{ - for (unsigned int d = 0; d < dim; ++d) - values (d) = d + 1.0; -} - -template -void test_boundary_values (const FiniteElement &fe, ConstraintMatrix &constraints) -{ - Triangulation triangulation; - GridGenerator::subdivided_hyper_cube (triangulation, 2); - DoFHandler dof_handler (triangulation); - dof_handler.distribute_dofs (fe); - BoundaryFunction boundary_function; - constraints.clear (); - VectorTools::project_boundary_values_curl_conforming (dof_handler, 0, boundary_function, 0, constraints); - constraints.close (); -} - -template -void test(unsigned order) -{ - deallog << "dim:" << dim << " order:" << order << "\t"; - - ConstraintMatrix constraints_fe, constraints_fes; - - { - FE_Nedelec<3> fe (order); - test_boundary_values (fe, constraints_fe); - } - - { - FESystem<3> fe(FE_Nedelec<3>(order),1); - test_boundary_values (fe, constraints_fes); - } - - if (constraints_fes.n_constraints() == constraints_fe.n_constraints()) - { - const IndexSet &lines = constraints_fes.get_local_lines (); - - for (unsigned i = 0; i < lines.n_elements(); ++i) - { - if (!constraints_fe.is_constrained(lines.nth_index_in_set(i))) - { - deallog << "Failed" << std::endl; - return; - } - - const std::vector > &c1 - = *constraints_fes.get_constraint_entries(lines.nth_index_in_set(i)); - const std::vector > &c2 - = *constraints_fe.get_constraint_entries(lines.nth_index_in_set(i)); - - for (size_t j = 0; j < c1.size(); ++j) - if ((c1[j].first != c2[j].first) || (fabs(c1[j].second - c2[j].second) > 1e-14)) - { - deallog << "Failed" << std::endl; - return; - } - } - } - else - { - deallog << "Failed" << std::endl; - return; - } - - deallog << "OK" << std::endl; -} - -int main () -{ - deallog << std::setprecision (2); - deallog.attach (logfile); - - test<2>(0); - test<2>(1); - test<2>(2); - - test<3>(0); - test<3>(1); - //test<3>(2); -} +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 2015 by the Alexander Grayver & deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// The test checks that project_boundary_values_curl_conforming +// works correctly for high-order FE_Nedelec elements used via +// FESystem. This requires the produced constraints to be the same +// for FE_Nedelec and FESystem(FE_Nedelec, 1). + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + +std::ofstream logfile ("output"); + +template +class BoundaryFunction: public Function +{ +public: + BoundaryFunction (); + virtual void vector_value (const Point &p, Vector &values) const; +}; + +template +BoundaryFunction::BoundaryFunction (): Function (dim) +{ +} + +template +void BoundaryFunction::vector_value (const Point &, + Vector &values) const +{ + for (unsigned int d = 0; d < dim; ++d) + values (d) = d + 1.0; +} + +template +void test_boundary_values (const FiniteElement &fe, ConstraintMatrix &constraints) +{ + Triangulation triangulation; + GridGenerator::subdivided_hyper_cube (triangulation, 2); + DoFHandler dof_handler (triangulation); + dof_handler.distribute_dofs (fe); + BoundaryFunction boundary_function; + constraints.clear (); + VectorTools::project_boundary_values_curl_conforming (dof_handler, 0, boundary_function, 0, constraints); + constraints.close (); +} + +template +void test(unsigned order) +{ + deallog << "dim:" << dim << " order:" << order << "\t"; + + ConstraintMatrix constraints_fe, constraints_fes; + + { + FE_Nedelec<3> fe (order); + test_boundary_values (fe, constraints_fe); + } + + { + FESystem<3> fe(FE_Nedelec<3>(order),1); + test_boundary_values (fe, constraints_fes); + } + + if (constraints_fes.n_constraints() == constraints_fe.n_constraints()) + { + const IndexSet &lines = constraints_fes.get_local_lines (); + + for (unsigned i = 0; i < lines.n_elements(); ++i) + { + if (!constraints_fe.is_constrained(lines.nth_index_in_set(i))) + { + deallog << "Failed" << std::endl; + return; + } + + const std::vector > &c1 + = *constraints_fes.get_constraint_entries(lines.nth_index_in_set(i)); + const std::vector > &c2 + = *constraints_fe.get_constraint_entries(lines.nth_index_in_set(i)); + + for (size_t j = 0; j < c1.size(); ++j) + if ((c1[j].first != c2[j].first) || (fabs(c1[j].second - c2[j].second) > 1e-14)) + { + deallog << "Failed" << std::endl; + return; + } + } + } + else + { + deallog << "Failed" << std::endl; + return; + } + + deallog << "OK" << std::endl; +} + +int main () +{ + deallog << std::setprecision (2); + deallog.attach (logfile); + + test<2>(0); + test<2>(1); + test<2>(2); + + test<3>(0); + test<3>(1); + //test<3>(2); +} diff --git a/tests/quick_tests/step-metis.cc b/tests/quick_tests/step-metis.cc index a0e31ab5c3..7bf9cb9e26 100644 --- a/tests/quick_tests/step-metis.cc +++ b/tests/quick_tests/step-metis.cc @@ -1,79 +1,79 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2013 - 2016 by the deal.II authors - * - * This file is part of the deal.II library. - * - * The deal.II library is free software; you can use it, redistribute - * it, and/or modify it under the terms of the GNU Lesser General - * Public License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE at - * the top level of the deal.II distribution. - * - * --------------------------------------------------------------------- - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -using namespace dealii; - -static const unsigned int dim = 2; - -int main (int argc, char **argv) -{ - try - { - Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); - { - parallel::shared::Triangulation triangulation(MPI_COMM_WORLD); - FE_Q fe(1); - DoFHandler dof_handler (triangulation); - - GridGenerator::hyper_cube (triangulation, -1, 1); - triangulation.refine_global (2); - dof_handler.distribute_dofs (fe); - IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); - deallog << locally_owned_dofs.n_elements() << std::endl; - dof_handler.clear (); - deallog << "OK" << std::endl; - } - } - - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} +/* --------------------------------------------------------------------- + * + * Copyright (C) 2013 - 2016 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +static const unsigned int dim = 2; + +int main (int argc, char **argv) +{ + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + parallel::shared::Triangulation triangulation(MPI_COMM_WORLD); + FE_Q fe(1); + DoFHandler dof_handler (triangulation); + + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (2); + dof_handler.distribute_dofs (fe); + IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); + deallog << locally_owned_dofs.n_elements() << std::endl; + dof_handler.clear (); + deallog << "OK" << std::endl; + } + } + + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/quick_tests/step-petsc.cc b/tests/quick_tests/step-petsc.cc index 7ffec9200f..419f0574ba 100644 --- a/tests/quick_tests/step-petsc.cc +++ b/tests/quick_tests/step-petsc.cc @@ -1,218 +1,218 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2013 - 2016 by the deal.II authors - * - * This file is part of the deal.II library. - * - * The deal.II library is free software; you can use it, redistribute - * it, and/or modify it under the terms of the GNU Lesser General - * Public License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE at - * the top level of the deal.II distribution. - * - * --------------------------------------------------------------------- - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -using namespace dealii; - -// Test that deal.II is working with PETSc by solving the Laplace's -// problem in 2d. -class LaplaceProblem -{ -public: - LaplaceProblem (); - void run (); - -private: - void setup_system (); - void assemble_system (); - void solve (); - - Triangulation<2> triangulation; - FE_Q<2> fe; - DoFHandler<2> dof_handler; - - PETScWrappers::SparseMatrix A; - PETScWrappers::MPI::Vector b, x; - ConstraintMatrix constraints; - - TableHandler output_table; -}; - -LaplaceProblem::LaplaceProblem () - : - fe (1), - dof_handler (triangulation) -{} - -void LaplaceProblem::setup_system () -{ - dof_handler.distribute_dofs (fe); - - constraints.clear (); - DoFTools::make_zero_boundary_constraints (dof_handler, constraints); - constraints.close (); - - A.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), - dof_handler.max_couplings_between_dofs()); - b.reinit (MPI_COMM_WORLD, dof_handler.n_dofs(), dof_handler.n_dofs()); - x.reinit (MPI_COMM_WORLD, dof_handler.n_dofs(), dof_handler.n_dofs()); - - // some output - output_table.add_value ("cells", triangulation.n_active_cells()); - output_table.add_value ("dofs", dof_handler.n_dofs()); -} - -void LaplaceProblem::assemble_system () -{ - QGauss<2> quadrature_formula(2); - - FEValues<2> fe_values (fe, quadrature_formula, - update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_A (dofs_per_cell, dofs_per_cell); - Vector cell_b (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - DoFHandler<2>::active_cell_iterator - cell = dof_handler.begin_active (), - endc = dof_handler.end (); - - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - cell_A = 0; - cell_b = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - constraints.distribute_local_to_global (cell_A, local_dof_indices, A); - constraints.distribute_local_to_global (cell_b, local_dof_indices, b); - } - - A.compress (VectorOperation::add); - b.compress (VectorOperation::add); -} - -void LaplaceProblem::solve () -{ - SolverControl solver_control (1e03, 1e-03); - PETScWrappers::SolverCG cg_solver (solver_control); - PETScWrappers::PreconditionBlockJacobi preconditioner (A); - cg_solver.solve (A, x, b, preconditioner); - - PETScWrappers::MPI::Vector res(x); - A.residual(res,x,b); - AssertThrow(res.l2_norm()<1e-3, - ExcInternalError()); -} - -void LaplaceProblem::run () -{ - GridGenerator::hyper_cube (triangulation, -1, 1); - - for (unsigned int c=0; c<5; ++c) - { - triangulation.refine_global (1); - setup_system (); - assemble_system (); - solve (); - } - - // finialise output - output_table.write_text (std::cout); - deallog << std::endl; -} - - -int main (int argc, char **argv) -{ - try - { - Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); - { - LaplaceProblem problem; - problem.run (); - deallog << "OK" << std::endl; - } - } - - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} +/* --------------------------------------------------------------------- + * + * Copyright (C) 2013 - 2016 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +// Test that deal.II is working with PETSc by solving the Laplace's +// problem in 2d. +class LaplaceProblem +{ +public: + LaplaceProblem (); + void run (); + +private: + void setup_system (); + void assemble_system (); + void solve (); + + Triangulation<2> triangulation; + FE_Q<2> fe; + DoFHandler<2> dof_handler; + + PETScWrappers::SparseMatrix A; + PETScWrappers::MPI::Vector b, x; + ConstraintMatrix constraints; + + TableHandler output_table; +}; + +LaplaceProblem::LaplaceProblem () + : + fe (1), + dof_handler (triangulation) +{} + +void LaplaceProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + constraints.clear (); + DoFTools::make_zero_boundary_constraints (dof_handler, constraints); + constraints.close (); + + A.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + b.reinit (MPI_COMM_WORLD, dof_handler.n_dofs(), dof_handler.n_dofs()); + x.reinit (MPI_COMM_WORLD, dof_handler.n_dofs(), dof_handler.n_dofs()); + + // some output + output_table.add_value ("cells", triangulation.n_active_cells()); + output_table.add_value ("dofs", dof_handler.n_dofs()); +} + +void LaplaceProblem::assemble_system () +{ + QGauss<2> quadrature_formula(2); + + FEValues<2> fe_values (fe, quadrature_formula, + update_values | + update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_A (dofs_per_cell, dofs_per_cell); + Vector cell_b (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + DoFHandler<2>::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_A = 0; + cell_b = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints.distribute_local_to_global (cell_A, local_dof_indices, A); + constraints.distribute_local_to_global (cell_b, local_dof_indices, b); + } + + A.compress (VectorOperation::add); + b.compress (VectorOperation::add); +} + +void LaplaceProblem::solve () +{ + SolverControl solver_control (1e03, 1e-03); + PETScWrappers::SolverCG cg_solver (solver_control); + PETScWrappers::PreconditionBlockJacobi preconditioner (A); + cg_solver.solve (A, x, b, preconditioner); + + PETScWrappers::MPI::Vector res(x); + A.residual(res,x,b); + AssertThrow(res.l2_norm()<1e-3, + ExcInternalError()); +} + +void LaplaceProblem::run () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + + for (unsigned int c=0; c<5; ++c) + { + triangulation.refine_global (1); + setup_system (); + assemble_system (); + solve (); + } + + // finialise output + output_table.write_text (std::cout); + deallog << std::endl; +} + + +int main (int argc, char **argv) +{ + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + LaplaceProblem problem; + problem.run (); + deallog << "OK" << std::endl; + } + } + + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/quick_tests/step-slepc.cc b/tests/quick_tests/step-slepc.cc index 8aab2313fa..b38d647f3c 100644 --- a/tests/quick_tests/step-slepc.cc +++ b/tests/quick_tests/step-slepc.cc @@ -1,261 +1,261 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2013 - 2016 by the deal.II authors - * - * This file is part of the deal.II library. - * - * The deal.II library is free software; you can use it, redistribute - * it, and/or modify it under the terms of the GNU Lesser General - * Public License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE at - * the top level of the deal.II distribution. - * - * --------------------------------------------------------------------- - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -using namespace dealii; - -// Test that deal.II is working with SLEPc by solving the Laplace's -// eigenspectrum problem in 2d. -class LaplaceEigenspectrumProblem -{ -public: - LaplaceEigenspectrumProblem (); - void run (); - -private: - void setup_system (); - void assemble_system (); - void solve (); - - Triangulation<2> triangulation; - FE_Q<2> fe; - DoFHandler<2> dof_handler; - - PETScWrappers::SparseMatrix A, B; - std::vector x; - std::vector lambda; - ConstraintMatrix constraints; - - TableHandler output_table; -}; - -LaplaceEigenspectrumProblem::LaplaceEigenspectrumProblem () - : - fe (1), - dof_handler (triangulation) -{} - -void LaplaceEigenspectrumProblem::setup_system () -{ - dof_handler.distribute_dofs (fe); - - constraints.clear (); - DoFTools::make_zero_boundary_constraints (dof_handler, constraints); - constraints.close (); - - A.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), - dof_handler.max_couplings_between_dofs()); - B.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), - dof_handler.max_couplings_between_dofs()); - - x.resize (1); - x[0].reinit (MPI_COMM_WORLD, dof_handler.n_dofs(), dof_handler.n_dofs()); - lambda.resize (1); - lambda[0] = 0.; - - // some output - output_table.add_value ("cells", triangulation.n_active_cells()); - output_table.add_value ("dofs", dof_handler.n_dofs()); -} - -void LaplaceEigenspectrumProblem::assemble_system () -{ - QGauss<2> quadrature_formula(2); - - FEValues<2> fe_values (fe, quadrature_formula, - update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_A (dofs_per_cell, dofs_per_cell); - FullMatrix cell_B (dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler<2>::active_cell_iterator - cell = dof_handler.begin_active (), - endc = dof_handler.end (); - - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - cell_A = 0; - cell_B = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - constraints.distribute_local_to_global (cell_A, local_dof_indices, A); - constraints.distribute_local_to_global (cell_B, local_dof_indices, B); - } - - A.compress (VectorOperation::add); - B.compress (VectorOperation::add); -} - -void LaplaceEigenspectrumProblem::solve () -{ - SolverControl solver_control (1000, 1e-10); - SLEPcWrappers::SolverArnoldi eigensolver (solver_control); - eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); - eigensolver.solve (A, B, lambda, x, x.size()); - - { - const double precision = 1e-7; - PETScWrappers::MPI::Vector Ax(x[0]), Bx(x[0]); - for (unsigned int i=0; i < x.size(); ++i) - { - B.vmult(Bx,x[i]); - - for (unsigned int j=0; j < x.size(); j++) - if (j!=i) - Assert( std::abs( x[j] * Bx )< precision, - ExcMessage("Eigenvectors " + - Utilities::int_to_string(i) + - " and " + - Utilities::int_to_string(j) + - " are not orthogonal!")); - - A.vmult(Ax,x[i]); - Ax.add(-1.0*lambda[i],Bx); - Assert (Ax.l2_norm() < precision, - ExcMessage("Returned vector " + - Utilities::int_to_string(i) + - " is not an eigenvector!")); - } - } - - - // some output - output_table.add_value ("lambda", lambda[0]); - output_table.add_value ("error", std::fabs(2.-lambda[0])); -} - -void LaplaceEigenspectrumProblem::run () -{ - const double radius = dealii::numbers::PI/2.; - GridGenerator::hyper_cube (triangulation, -radius, radius); - - // set the old eigenvalue to a silly number. - double old_lambda = 1000; - - for (unsigned int c=0; c<5; ++c) - { - // obtain numerical result - triangulation.refine_global (1); - setup_system (); - assemble_system (); - solve (); - - // check energy convergence with previous result - AssertThrow (lambda[0] +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +// Test that deal.II is working with SLEPc by solving the Laplace's +// eigenspectrum problem in 2d. +class LaplaceEigenspectrumProblem +{ +public: + LaplaceEigenspectrumProblem (); + void run (); + +private: + void setup_system (); + void assemble_system (); + void solve (); + + Triangulation<2> triangulation; + FE_Q<2> fe; + DoFHandler<2> dof_handler; + + PETScWrappers::SparseMatrix A, B; + std::vector x; + std::vector lambda; + ConstraintMatrix constraints; + + TableHandler output_table; +}; + +LaplaceEigenspectrumProblem::LaplaceEigenspectrumProblem () + : + fe (1), + dof_handler (triangulation) +{} + +void LaplaceEigenspectrumProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + constraints.clear (); + DoFTools::make_zero_boundary_constraints (dof_handler, constraints); + constraints.close (); + + A.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + B.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + + x.resize (1); + x[0].reinit (MPI_COMM_WORLD, dof_handler.n_dofs(), dof_handler.n_dofs()); + lambda.resize (1); + lambda[0] = 0.; + + // some output + output_table.add_value ("cells", triangulation.n_active_cells()); + output_table.add_value ("dofs", dof_handler.n_dofs()); +} + +void LaplaceEigenspectrumProblem::assemble_system () +{ + QGauss<2> quadrature_formula(2); + + FEValues<2> fe_values (fe, quadrature_formula, + update_values | + update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_A (dofs_per_cell, dofs_per_cell); + FullMatrix cell_B (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler<2>::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_A = 0; + cell_B = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints.distribute_local_to_global (cell_A, local_dof_indices, A); + constraints.distribute_local_to_global (cell_B, local_dof_indices, B); + } + + A.compress (VectorOperation::add); + B.compress (VectorOperation::add); +} + +void LaplaceEigenspectrumProblem::solve () +{ + SolverControl solver_control (1000, 1e-10); + SLEPcWrappers::SolverArnoldi eigensolver (solver_control); + eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); + eigensolver.solve (A, B, lambda, x, x.size()); + + { + const double precision = 1e-7; + PETScWrappers::MPI::Vector Ax(x[0]), Bx(x[0]); + for (unsigned int i=0; i < x.size(); ++i) + { + B.vmult(Bx,x[i]); + + for (unsigned int j=0; j < x.size(); j++) + if (j!=i) + Assert( std::abs( x[j] * Bx )< precision, + ExcMessage("Eigenvectors " + + Utilities::int_to_string(i) + + " and " + + Utilities::int_to_string(j) + + " are not orthogonal!")); + + A.vmult(Ax,x[i]); + Ax.add(-1.0*lambda[i],Bx); + Assert (Ax.l2_norm() < precision, + ExcMessage("Returned vector " + + Utilities::int_to_string(i) + + " is not an eigenvector!")); + } + } + + + // some output + output_table.add_value ("lambda", lambda[0]); + output_table.add_value ("error", std::fabs(2.-lambda[0])); +} + +void LaplaceEigenspectrumProblem::run () +{ + const double radius = dealii::numbers::PI/2.; + GridGenerator::hyper_cube (triangulation, -radius, radius); + + // set the old eigenvalue to a silly number. + double old_lambda = 1000; + + for (unsigned int c=0; c<5; ++c) + { + // obtain numerical result + triangulation.refine_global (1); + setup_system (); + assemble_system (); + solve (); + + // check energy convergence with previous result + AssertThrow (lambda[0] -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -using namespace dealii; - -// Test that deal.II is working with Trilinos by solving the Laplace's -// problem in 2d. -class LaplaceProblem -{ -public: - LaplaceProblem (); - void run (); - -private: - void setup_system (); - void assemble_system (); - void solve (); - - Triangulation<2> triangulation; - FE_Q<2> fe; - DoFHandler<2> dof_handler; - - TrilinosWrappers::SparseMatrix A; - TrilinosWrappers::MPI::Vector b, x; - ConstraintMatrix constraints; - - TableHandler output_table; -}; - -LaplaceProblem::LaplaceProblem () - : - fe (1), - dof_handler (triangulation) -{} - -void LaplaceProblem::setup_system () -{ - dof_handler.distribute_dofs (fe); - - constraints.clear (); - DoFTools::make_zero_boundary_constraints (dof_handler, constraints); - constraints.close (); - - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false); - - A.reinit (dsp); - b.reinit (complete_index_set(dof_handler.n_dofs())); - x.reinit (complete_index_set(dof_handler.n_dofs())); - - // some output - output_table.add_value ("cells", triangulation.n_active_cells()); - output_table.add_value ("dofs", dof_handler.n_dofs()); -} - -void LaplaceProblem::assemble_system () -{ - QGauss<2> quadrature_formula(2); - - FEValues<2> fe_values (fe, quadrature_formula, - update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_A (dofs_per_cell, dofs_per_cell); - Vector cell_b (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - DoFHandler<2>::active_cell_iterator - cell = dof_handler.begin_active (), - endc = dof_handler.end (); - - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - cell_A = 0; - cell_b = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - constraints.distribute_local_to_global (cell_A, local_dof_indices, A); - constraints.distribute_local_to_global (cell_b, local_dof_indices, b); - } - - A.compress (VectorOperation::add); - b.compress (VectorOperation::add); -} - -void LaplaceProblem::solve () -{ - SolverControl solver_control (1e03, 1e-03); - TrilinosWrappers::SolverCG cg_solver (solver_control); - TrilinosWrappers::PreconditionBlockJacobi preconditioner; - preconditioner.initialize (A); - cg_solver.solve (A, x, b, preconditioner); - - TrilinosWrappers::MPI::Vector res(x); - A.residual(res,x,b); - AssertThrow(res.l2_norm()<1e-3, - ExcInternalError()); -} - -void LaplaceProblem::run () -{ - GridGenerator::hyper_cube (triangulation, -1, 1); - - for (unsigned int c=0; c<5; ++c) - { - triangulation.refine_global (1); - setup_system (); - assemble_system (); - solve (); - } - - // finialise output - output_table.write_text (std::cout); - deallog << std::endl; -} - - -int main (int argc, char **argv) -{ - try - { - Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); - { - LaplaceProblem problem; - problem.run (); - deallog << "OK" << std::endl; - } - } - - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} +/* --------------------------------------------------------------------- + * + * Copyright (C) 2013 - 2016 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +// Test that deal.II is working with Trilinos by solving the Laplace's +// problem in 2d. +class LaplaceProblem +{ +public: + LaplaceProblem (); + void run (); + +private: + void setup_system (); + void assemble_system (); + void solve (); + + Triangulation<2> triangulation; + FE_Q<2> fe; + DoFHandler<2> dof_handler; + + TrilinosWrappers::SparseMatrix A; + TrilinosWrappers::MPI::Vector b, x; + ConstraintMatrix constraints; + + TableHandler output_table; +}; + +LaplaceProblem::LaplaceProblem () + : + fe (1), + dof_handler (triangulation) +{} + +void LaplaceProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + constraints.clear (); + DoFTools::make_zero_boundary_constraints (dof_handler, constraints); + constraints.close (); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false); + + A.reinit (dsp); + b.reinit (complete_index_set(dof_handler.n_dofs())); + x.reinit (complete_index_set(dof_handler.n_dofs())); + + // some output + output_table.add_value ("cells", triangulation.n_active_cells()); + output_table.add_value ("dofs", dof_handler.n_dofs()); +} + +void LaplaceProblem::assemble_system () +{ + QGauss<2> quadrature_formula(2); + + FEValues<2> fe_values (fe, quadrature_formula, + update_values | + update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_A (dofs_per_cell, dofs_per_cell); + Vector cell_b (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + DoFHandler<2>::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_A = 0; + cell_b = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints.distribute_local_to_global (cell_A, local_dof_indices, A); + constraints.distribute_local_to_global (cell_b, local_dof_indices, b); + } + + A.compress (VectorOperation::add); + b.compress (VectorOperation::add); +} + +void LaplaceProblem::solve () +{ + SolverControl solver_control (1e03, 1e-03); + TrilinosWrappers::SolverCG cg_solver (solver_control); + TrilinosWrappers::PreconditionBlockJacobi preconditioner; + preconditioner.initialize (A); + cg_solver.solve (A, x, b, preconditioner); + + TrilinosWrappers::MPI::Vector res(x); + A.residual(res,x,b); + AssertThrow(res.l2_norm()<1e-3, + ExcInternalError()); +} + +void LaplaceProblem::run () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + + for (unsigned int c=0; c<5; ++c) + { + triangulation.refine_global (1); + setup_system (); + assemble_system (); + solve (); + } + + // finialise output + output_table.write_text (std::cout); + deallog << std::endl; +} + + +int main (int argc, char **argv) +{ + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + LaplaceProblem problem; + problem.run (); + deallog << "OK" << std::endl; + } + } + + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +}