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