]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Consistently use UNIX line endings. 5425/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Wed, 8 Nov 2017 17:42:04 +0000 (10:42 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 8 Nov 2017 17:42:04 +0000 (10:42 -0700)
12 files changed:
tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc
tests/adolc/basic_functions_jacobian_taped.cc
tests/adolc/basic_functions_jacobian_tapeless.cc
tests/distributed_grids/intergrid_transfer_representation_parallel.cc
tests/fe/fe_project_2d.cc
tests/fe/fe_project_3d.cc
tests/numerics/no_flux_13.cc
tests/numerics/project_bv_curl_conf_02.cc
tests/quick_tests/step-metis.cc
tests/quick_tests/step-petsc.cc
tests/quick_tests/step-slepc.cc
tests/quick_tests/step-trilinos.cc

index 9b3b2cb0328d0a176deef03b74714e19929678e7..e5ae2c5ee8ab71bce01e57962680c7f610c6f2cb 100644 (file)
-// ---------------------------------------------------------------------\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
+}
index 76a4c3b261a0bd8faa524e8df07f578432319a78..740931025d39add9bec8e99dcfc8cec2e1e1492c 100644 (file)
-//\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;
+}
index 8efccb0e0d751a0c7f0520df87ea6110ac6fd5b4..ab06cdf950570c15da30e8f7605fc9209c82f31d 100644 (file)
@@ -1,82 +1,82 @@
-// ---------------------------------------------------------------------\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;
+}
index cbb3adc4d9c0009f8110e8588e86201dc7668a7f..393d87ff0dd6526ac1b94da7022122f0be5ca821 100644 (file)
-// ---------------------------------------------------------------------\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);
+    }
+}
index 7d0b15c02188bcc0ad6bab9ce89d2e30ed9dc5ca..696b68a92337611aca77bc70091ac29375c599ce 100644 (file)
-//\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;
+}
index 590b25d486064b76c1c916ba07d0ce58b7b367f9..864b212be2345bfdbd56b8c4ec3f65b1fadce50d 100644 (file)
-//\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;
+}
index 0f472b4dbafa5e932d41ec0618fe7e3b4fba1351..6e32dcdcb75080638ed996737d6559ca27cc8502 100644 (file)
@@ -1,71 +1,71 @@
-// ---------------------------------------------------------------------\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> ();
+}
index 118514f036c15cf690c138da5a36df92c3d82045..46552e0737ecda8acfa2e755d5eaf4f24494add6 100644 (file)
-// ---------------------------------------------------------------------\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);
+}
index a0e31ab5c388d94fa51c2bab0de4e71947b1a15f..7bf9cb9e268ffc72702503bd903d60e5fbc35d3d 100644 (file)
@@ -1,79 +1,79 @@
-/* ---------------------------------------------------------------------\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;
+}
index 7ffec9200f50863fcbab176c3a97b96721d41e62..419f0574ba8043cc85d285b1f392be07e979530b 100644 (file)
-/* ---------------------------------------------------------------------\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;
+}
index 8aab2313fa161d60dba370cfbb5bc03cb2e346f6..b38d647f3cbac8fa56f7eb4275b53e1250adc733 100644 (file)
-/* ---------------------------------------------------------------------\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;
+}
index 29a1c6235fd238bca0776bd0956e157d12ab4f33..cc4a0dd69d03774426a4fef30453e13a160364d0 100644 (file)
-/* ---------------------------------------------------------------------\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;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.