]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Compute restriction and prolongation matrices for DGQ elements also if dim<spacedim.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 9 Feb 2011 15:30:57 +0000 (15:30 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 9 Feb 2011 15:30:57 +0000 (15:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@23314 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/source/fe/fe_dgq.cc
tests/codim_one/fe_dgq_prolongation_01.cc [new file with mode: 0644]
tests/codim_one/fe_dgq_prolongation_01/cmp/generic [new file with mode: 0644]
tests/codim_one/solution_transfer_01.cc [new file with mode: 0644]
tests/codim_one/solution_transfer_01/cmp/generic [new file with mode: 0644]

index b8b1095e10014050379eb73b5bab7dc2b522b448..f52355a0ee3fc7b9d73d4db27b5ca3ceb3dbce3e 100644 (file)
@@ -78,6 +78,15 @@ should be fixed now.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: Prolongation and restriction matrices were not computed at all
+for elements of type FE_DGQ if <code>dim@<spacedim</code>. Consequently,
+consumers of this information, such as the SolutionTransfer class or
+the DoFCellAccess::set_dof_values_by_interpolation function did not
+work either and simply returned zero results. This is now fixed.
+<br>
+(M. Sebastian Pauletti, Wolfgang Bangerth, 2011/02/09)
+</li>
+
 <li> Fixed: When refining a mesh with codimension one, edges were refined using
 the same manifold description as adjacent cells, but this ignored that a
 boundary indicator might have been purposefully set for edges that are truly at
index 3577c85c5615284544f193e9a14c7f6e636bcb3c..11d3216cca11ad3c85d968f7d17054c015871d2a 100644 (file)
@@ -142,12 +142,32 @@ FE_DGQ<dim, spacedim>::FE_DGQ (const unsigned int degree)
                                   // restriction and prolongation
                                   // matrices to the right sizes
   this->reinit_restriction_and_prolongation_matrices();
-                                  // Fill prolongation matrices with embedding operators
-  if (dim ==spacedim)
-  FETools::compute_embedding_matrices (*this, this->prolongation);
-                                  // Fill restriction matrices with L2-projection
-  if (dim ==spacedim)
-  FETools::compute_projection_matrices (*this, this->restriction);
+
+                                  // Fill prolongation matrices with
+                                  // embedding operators and
+                                  // restriction with L2 projection
+                                  //
+                                  // we can call the respective
+                                  // functions in the case
+                                  // dim==spacedim, but they are not
+                                  // currently implemented for the
+                                  // codim-1 case. there's no harm
+                                  // here, though, since these
+                                  // matrices are independent of the
+                                  // space dimension and so we can
+                                  // copy them from the respective
+                                  // elements (not cheap, but works)
+  if (dim == spacedim)
+    {
+      FETools::compute_embedding_matrices (*this, this->prolongation);
+      FETools::compute_projection_matrices (*this, this->restriction);
+    }
+  else
+    {
+      FE_DGQ<dim> tmp (degree);
+      this->prolongation = tmp.prolongation;
+      this->restriction  = tmp.restriction;
+    }
 
 
                                   // finally fill in support points
diff --git a/tests/codim_one/fe_dgq_prolongation_01.cc b/tests/codim_one/fe_dgq_prolongation_01.cc
new file mode 100644 (file)
index 0000000..567cdf9
--- /dev/null
@@ -0,0 +1,45 @@
+//----------------------------  fe_dgq_prolongation_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  fe_dgq_prolongation_01.cc  ---------------------------
+
+// this is the root cause for solution_tranfer_01: the prolongation
+// matrices for FE_DGQ<dim-1,dim> were not computed at all
+
+#include "../tests.h"
+#include <fe/fe_dgq.h>
+
+
+#include <fstream>
+
+using namespace dealii;
+
+
+int main ()
+{
+  std::ofstream logfile("fe_dgq_prolongation_01/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  const unsigned int spacedim = 2;
+  const unsigned int dim = spacedim-1;
+
+  for (unsigned int degree=0; degree<3; ++degree)
+    {
+      deallog << "Degree=" << degree << std::endl;
+      FE_DGQ<dim,spacedim> fe(degree);
+      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+       for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+         deallog << fe.get_prolongation_matrix(0)(i,j) << std::endl;
+    }
+
+  return 0;
+}
diff --git a/tests/codim_one/fe_dgq_prolongation_01/cmp/generic b/tests/codim_one/fe_dgq_prolongation_01/cmp/generic
new file mode 100644 (file)
index 0000000..05d8ce9
--- /dev/null
@@ -0,0 +1,18 @@
+
+DEAL::Degree=0
+DEAL::1.00000
+DEAL::Degree=1
+DEAL::1.00000
+DEAL::0
+DEAL::0.500000
+DEAL::0.500000
+DEAL::Degree=2
+DEAL::1.00000
+DEAL::0
+DEAL::0
+DEAL::0.375000
+DEAL::0.750000
+DEAL::-0.125000
+DEAL::0
+DEAL::1.00000
+DEAL::0
diff --git a/tests/codim_one/solution_transfer_01.cc b/tests/codim_one/solution_transfer_01.cc
new file mode 100644 (file)
index 0000000..944ad7a
--- /dev/null
@@ -0,0 +1,100 @@
+//----------------------------  solution_transfer_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  solution_transfer_01.cc  ---------------------------
+
+
+/*
+  Bug for solution transfer using DG in surfaces
+
+
+  set solution = 1
+  refine
+  do solution transfer
+  And the magically, solution becomes 0
+
+  This surprisingly turned out to be a problem with the prolongation
+  matrices of DGQ in codim-1. (See the fe_dgq_prolongation_01 test.)
+*/
+
+#include "../tests.h"
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <grid/grid_tools.h>
+#include <grid/grid_refinement.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_dgq.h>
+
+#include <numerics/vectors.h>
+#include <numerics/solution_transfer.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+
+int main ()
+{
+  std::ofstream logfile("solution_transfer_01/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  const unsigned int spacedim = 2;
+  const unsigned int dim = spacedim-1;
+
+  Triangulation<dim,spacedim> boundary_mesh;
+
+                                  // create a line in 2d
+  GridGenerator::hyper_cube (boundary_mesh);
+
+                                  // attach a piecewise constant
+                                  // element to this one cell
+  FE_DGQ<dim,spacedim> fe(0);
+  DoFHandler<dim,spacedim>  dh (boundary_mesh);
+  dh.distribute_dofs(fe);
+
+  Vector<double> solution(dh.n_dofs());
+  solution = 1.0;
+
+  deallog << "Old values:" << std::endl;
+  for (unsigned int i=0; i<solution.size(); i++)
+    deallog << solution(i) << endl;
+
+
+                                  // Do some refinement
+  boundary_mesh.begin_active()->set_refine_flag ();
+
+  SolutionTransfer<dim, Vector<double>, DoFHandler<dim,spacedim> >
+    soltrans(dh);
+
+  boundary_mesh.prepare_coarsening_and_refinement();
+
+  soltrans.prepare_for_coarsening_and_refinement(solution);
+  boundary_mesh.execute_coarsening_and_refinement ();
+
+  dh.distribute_dofs(fe);
+
+                                  // get the interpolated solution
+                                  // back
+  Vector<double> tmp(dh.n_dofs());
+  tmp = 2;
+  soltrans.interpolate(solution, tmp);
+
+  deallog << "New values:" << std::endl;
+  for (unsigned int i=0; i<tmp.size(); i++)
+    deallog << tmp(i) << endl;
+
+  return 0;
+}
diff --git a/tests/codim_one/solution_transfer_01/cmp/generic b/tests/codim_one/solution_transfer_01/cmp/generic
new file mode 100644 (file)
index 0000000..8408606
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::Old values:
+DEAL::1.00000
+DEAL::New values:
+DEAL::1.00000
+DEAL::1.00000

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.