]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merged from trunk r32177 through r32196.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 12 Jan 2014 10:54:07 +0000 (10:54 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 12 Jan 2014 10:54:07 +0000 (10:54 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32197 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/examples/step-44/step-44.cc
deal.II/include/deal.II/fe/fe.h
deal.II/include/deal.II/meshworker/simple.h
deal.II/source/distributed/grid_refinement.cc
deal.II/source/fe/fe.cc
deal.II/source/fe/fe_poly_tensor.cc
deal.II/source/fe/fe_system.cc
deal.II/source/fe/fe_values.cc

index 8b705acd64638047a5739a5365ef347914720c57..4c21c7596b20894ca312770d5482e0e5ea32cd27 100644 (file)
@@ -1,7 +1,7 @@
  // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 2013 by the deal.II authors
+// Copyright (C) 2013, 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -85,6 +85,27 @@ inconvenience this causes.
 
 <ol>
 
+  <li> Fixed: There were several bugs in functions like
+  FEValues::get_function_values() where the code did not properly handle the
+  case of FE_Nothing. This is now fixed.
+  <br>
+  (Wolfgang Bangerth, 2014/01/08)
+  </li>
+
+  <li> Fixed: DataOut got confused in some situations where one uses FE_Nothing.
+  This is now fixed.
+  <br>
+  (Minh Do-Quang, Wolfgang Bangerth, 2014/01/08)
+  </li>
+
+  <li> Fixed: FESystem::get_interpolation_matrix, a function that is among
+  other places used by SolutionTransfer, had a bug that prevented it from
+  running correctly in some situations where one uses FE_Nothing.
+  This is now fixed.
+  <br>
+  (Minh Do-Quang, Wolfgang Bangerth, 2014/01/08)
+  </li>
+
   <li> Improved: When you call WorkStream::run with an empty function object
   for the copier, operations on individual cells are essentially all independent.
   In other words, you have a massively parallel collection of jobs. In this
index 5e2764781041e8cae4dd39ebab339fd71a7f383a..f7b59a2f17bc2a7f8070ea8cd1bde5e6d3b76cac 100644 (file)
@@ -15,7 +15,7 @@
  *
  * ---------------------------------------------------------------------
 
-/*
+ *
  * Authors: Jean-Paul Pelteret, University of Cape Town,
  *          Andrew McBride, University of Erlangen-Nuremberg, 2010
  */
index 40b4c742c5964f283369d647f51446665effc6c8..fda2a7fef32a94f6213cf3d0114f88ce4a905cac 100644 (file)
@@ -931,6 +931,19 @@ public:
   std::pair<unsigned int, unsigned int>
   face_system_to_component_index (const unsigned int index) const;
 
+  /**
+   * For faces with non-standard face_orientation in 3D, the dofs on faces
+   * (quads) have to be permuted in order to be combined with the correct
+   * shape functions. Given a local dof @p index on a quad, return the local
+   * index, if the face has non-standard face_orientation, face_flip or
+   * face_rotation. In 2D and 1D there is no need for permutation and
+   * consequently an exception is thrown.
+   */
+  unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index,
+                                                           const bool face_orientation,
+                                                           const bool face_flip,
+                                                           const bool face_rotation) const;
+
   /**
    * Given an index in the natural ordering of indices on a face, return the
    * index of the same degree of freedom on the cell.
@@ -989,20 +1002,7 @@ public:
                                    const bool face_orientation = true,
                                    const bool face_flip        = false,
                                    const bool face_rotation    = false) const;
-
-  /**
-   * For faces with non-standard face_orientation in 3D, the dofs on faces
-   * (quads) have to be permuted in order to be combined with the correct
-   * shape functions. Given a local dof @p index on a quad, return the local
-   * index, if the face has non-standard face_orientation, face_flip or
-   * face_rotation. In 2D and 1D there is no need for permutation and
-   * consequently an exception is thrown.
-   */
-  unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index,
-                                                           const bool face_orientation,
-                                                           const bool face_flip,
-                                                           const bool face_rotation) const;
-
   /**
    * For lines with non-standard line_orientation in 3D, the dofs on lines
    * have to be permuted in order to be combined with the correct shape
index 224729da485d29a235ddda1cf043a6efe7b5c08c..ed692c57f3db468243b0398f81d93f7ae8e3dd1a 100644 (file)
@@ -975,18 +975,13 @@ namespace MeshWorker
             if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
                 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
               {
-                if (mg_constrained_dofs->set_boundary_values())
-                  {
-                    if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                         !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
-                        ||
-                        (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                         mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
-                         i1[j] == i2[k]))
-                      G.add(i1[j], i2[k], M(j,k));
-                  }
-                else
-                  G.add(i1[j], i2[k], M(j,k));
+               if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                    !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+                   ||
+                   (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                    mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                    i1[j] == i2[k]))
+                 G.add(i1[j], i2[k], M(j,k));
               }
     }
 
@@ -1010,18 +1005,13 @@ namespace MeshWorker
             if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
                 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
               {
-                if (mg_constrained_dofs->set_boundary_values())
-                  {
-                    if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                         !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
-                        ||
-                        (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
-                         mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
-                         i1[j] == i2[k]))
-                      G.add(i1[j], i2[k], M(k,j));
-                  }
-                else
-                  G.add(i1[j], i2[k], M(k,j));
+               if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                    !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+                   ||
+                   (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                    mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                    i1[j] == i2[k]))
+                 G.add(i1[j], i2[k], M(k,j));
               }
     }
 
index 9b89144f91a705d6e0c99474bd44318b6ec543df..f7f185ff4e1289f53ce5bf783c16f2e625fdc44b 100644 (file)
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 2000 - 2013 by the deal.II authors
+// Copyright (C) 2000 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -315,28 +315,15 @@ namespace
           else
             interesting_range[0] = interesting_range[1] = test_threshold;
 
-          // terminate the iteration
-          // after 10 go-arounds. this
-          // is necessary because
-          // oftentimes error
-          // indicators on cells have
-          // exactly the same value,
-          // and so there may not be a
-          // particular value that cuts
-          // the indicators in such a
-          // way that we can achieve
-          // the desired number of
-          // cells. using a max of 10
-          // iterations means that we
-          // terminate the iteration
-          // after 10 steps if the
-          // indicators were perfectly
-          // badly distributed, and we
-          // make at most a mistake of
-          // 1/2^10 in the number of
-          // cells flagged if
-          // indicators are perfectly
-          // equidistributed
+          // terminate the iteration after 25 go-arounds. this is necessary
+          // because oftentimes error indicators on cells have exactly the
+          // same value, and so there may not be a particular value that cuts
+          // the indicators in such a way that we can achieve the desired
+          // number of cells. using a maximal number of iterations means that
+          // we terminate the iteration after a fixed number N of steps if the
+          // indicators were perfectly badly distributed, and we make at most
+          // a mistake of 1/2^N in the number of cells flagged if indicators
+          // are perfectly equidistributed
           ++iteration;
           if (iteration == 25)
             interesting_range[0] = interesting_range[1] = test_threshold;
@@ -400,8 +387,9 @@ namespace
   {
     /**
      * Compute a threshold value so
-     * that exactly n_target_cells have
-     * a value that is larger.
+     * that the error accumulated over all criteria[i] so that
+     *     criteria[i] > threshold
+     * is larger than target_error.
      */
     template <typename number>
     number
index f8c048226103d1d9ec315054e49b996456cf46e3..43db5378c1d334b6b3f11188268a3ef445166f27 100644 (file)
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -557,7 +557,7 @@ face_to_cell_index (const unsigned int face_index,
           ExcIndexRange(face_index, 0, this->dofs_per_face));
   Assert (face < GeometryInfo<dim>::faces_per_cell,
           ExcIndexRange(face, 0, GeometryInfo<dim>::faces_per_cell));
-
 //TODO: we could presumably solve the 3d case below using the
 // adjust_quad_dof_index_for_face_orientation_table field. for the
 // 2d case, we can't use adjust_line_dof_index_for_line_orientation_table
@@ -566,7 +566,7 @@ face_to_cell_index (const unsigned int face_index,
 // DoFTools::make_periodicity_constraints, for example). so we
 // would need to either fill this field, or rely on derived classes
 // implementing this function, as we currently do
-
   // see the function's documentation for an explanation of this
   // assertion -- in essence, derived classes have to implement
   // an overloaded version of this function if we are to use any
@@ -577,7 +577,7 @@ face_to_cell_index (const unsigned int face_index,
                         "Rather, the derived class you are using must provide "
                         "an overloaded version but apparently hasn't done so. See "
                         "the documentation of this function for more information."));
-
   // we need to distinguish between DoFs on vertices, lines and in 3d quads.
   // do so in a sequence of if-else statements
   if (face_index < this->first_face_line_index)
@@ -587,7 +587,7 @@ face_to_cell_index (const unsigned int face_index,
       // along with the number of the DoF on this vertex
       const unsigned int face_vertex         = face_index / this->dofs_per_vertex;
       const unsigned int dof_index_on_vertex = face_index % this->dofs_per_vertex;
-
       // then get the number of this vertex on the cell and translate
       // this to a DoF number on the cell
       return (GeometryInfo<dim>::face_to_cell_vertices(face, face_vertex,
@@ -604,10 +604,10 @@ face_to_cell_index (const unsigned int face_index,
       // do the same kind of translation as before. we need to only consider
       // DoFs on the lines, i.e., ignoring those on the vertices
       const unsigned int index = face_index - this->first_face_line_index;
-
       const unsigned int face_line         = index / this->dofs_per_line;
       const unsigned int dof_index_on_line = index % this->dofs_per_line;
-
       return (this->first_line_index
               + GeometryInfo<dim>::face_to_cell_lines(face, face_line,
                                                       face_orientation,
@@ -621,19 +621,19 @@ face_to_cell_index (const unsigned int face_index,
     // DoF is on a quad
     {
       Assert (dim >= 3, ExcInternalError());
-
       // ignore vertex and line dofs
       const unsigned int index = face_index - this->first_face_quad_index;
-
       return (this->first_quad_index
               + face * this->dofs_per_quad
               + index);
     }
 }
-
-
-
-
 template <int dim, int spacedim>
 unsigned int
 FiniteElement<dim,spacedim>::adjust_quad_dof_index_for_face_orientation (const unsigned int index,
@@ -1248,6 +1248,11 @@ FiniteElement<dim,spacedim>::compute_2nd (
           & update_hessians,
           ExcInternalError());
 
+  // there is nothing to do if there are no degrees of freedom (e.g., in an
+  // FE_Nothing)
+  if (this->dofs_per_cell == 0)
+    return;
+
 // make sure we have as many entries as there are nonzero components
 //  Assert (data.shape_hessians.size() ==
 //          std::accumulate (n_nonzero_components_table.begin(),
index 9f6460c363d50fa3cc67440aa019912fa7b86628..4735a210ef462b369c7accc784b3a731637de7e3 100644 (file)
@@ -585,9 +585,9 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
     }
 
   if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
-    this->compute_2nd (mapping, cell,
-                       typename QProjector<dim>::DataSetDescriptor().cell(),
-                       mapping_data, fe_data, data);
+       this->compute_2nd (mapping, cell,
+                          typename QProjector<dim>::DataSetDescriptor().cell(),
+                          mapping_data, fe_data, data);
 }
 
 
@@ -791,7 +791,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
     }
 
   if (flags & update_hessians)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+       this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
 }
 
 
@@ -989,7 +989,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 
               for (unsigned int k = 0; k < n_q_points; ++k)
                 for (unsigned int d = 0; d < dim; ++d)
-                  data.shape_gradients[first + d][k] = shape_grads1[k][d];
+                  data.shape_gradients[first + d][k] = sign_change[i] * shape_grads1[k][d];
 
               break;
             }
@@ -1001,7 +1001,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
     }
 
   if (flags & update_hessians)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+       this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
 }
 
 
index bf70ec76a66a4172825259726658b8a54a4dfe91..8683136e667d7a1234e99679f203581a622923a5 100644 (file)
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 1999 - 2013 by the deal.II authors
+// Copyright (C) 1999 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -528,10 +528,18 @@ FESystem<dim,spacedim>::get_interpolation_matrix (
   const FiniteElement<dim,spacedim> &x_source_fe,
   FullMatrix<double>           &interpolation_matrix) const
 {
-  Assert (interpolation_matrix.m() == this->dofs_per_cell,
+  // check that the size of the matrices is correct. for historical
+  // reasons, if you call matrix.reinit(8,0), it sets the sizes
+  // to m==n==0 internally. this may happen when we use a FE_Nothing,
+  // so write the test in a more lenient way
+  Assert ((interpolation_matrix.m() == this->dofs_per_cell)
+         ||
+         (x_source_fe.dofs_per_cell == 0),
           ExcDimensionMismatch (interpolation_matrix.m(),
                                 this->dofs_per_cell));
-  Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell,
+  Assert ((interpolation_matrix.n() == x_source_fe.dofs_per_cell)
+         ||
+         (this->dofs_per_cell == 0),
           ExcDimensionMismatch (interpolation_matrix.m(),
                                 x_source_fe.dofs_per_cell));
 
index 066c79c4c870fa4b8439a46357df6d2d1e2baf98..aaf9fe0a0a746d355bc9e6670484848e0bc5a565 100644 (file)
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -2242,9 +2242,18 @@ namespace internal
                       const bool quadrature_points_fastest  = false,
                       const unsigned int component_multiple = 1)
   {
+    // initialize with zero
+    for (unsigned int i=0; i<values.size(); ++i)
+      std::fill_n (values[i].begin(), values[i].size(),
+                   typename VectorType::value_type());
+
+    // see if there the current cell has DoFs at all, and if not
+    // then there is nothing else to do.
     const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
-                                             shape_values.n_cols() : 0;
+    if (dofs_per_cell == 0)
+      return;
+    
+    const unsigned int n_quadrature_points = shape_values.n_cols();
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2262,11 +2271,6 @@ namespace internal
           AssertDimension (values[i].size(), result_components);
       }
 
-    // initialize with zero
-    for (unsigned int i=0; i<values.size(); ++i)
-      std::fill_n (values[i].begin(), values[i].size(),
-                   typename VectorType::value_type());
-
     // add up contributions of trial functions.  now check whether the shape
     // function is primitive or not. if it is, then set its only non-zero
     // component, otherwise loop over components
@@ -2369,9 +2373,19 @@ namespace internal
                            const bool quadrature_points_fastest  = false,
                            const unsigned int component_multiple = 1)
   {
+    // initialize with zero
+    for (unsigned int i=0; i<derivatives.size(); ++i)
+      std::fill_n (derivatives[i].begin(), derivatives[i].size(),
+                   Tensor<order,spacedim>());
+
+    // see if there the current cell has DoFs at all, and if not
+    // then there is nothing else to do.
     const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
-                                             shape_derivatives[0].size() : 0;
+    if (dofs_per_cell == 0)
+      return;
+
+
+    const unsigned int n_quadrature_points = shape_derivatives[0].size();
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2389,11 +2403,6 @@ namespace internal
           AssertDimension (derivatives[i].size(), result_components);
       }
 
-    // initialize with zero
-    for (unsigned int i=0; i<derivatives.size(); ++i)
-      std::fill_n (derivatives[i].begin(), derivatives[i].size(),
-                   Tensor<order,spacedim>());
-
     // add up contributions of trial functions.  now check whether the shape
     // function is primitive or not. if it is, then set its only non-zero
     // component, otherwise loop over components
@@ -2485,9 +2494,19 @@ namespace internal
                           const bool quadrature_points_fastest  = false,
                           const unsigned int component_multiple = 1)
   {
+    // initialize with zero
+    for (unsigned int i=0; i<laplacians.size(); ++i)
+      std::fill_n (laplacians[i].begin(), laplacians[i].size(),
+                   typename VectorType::value_type());
+
+    // see if there the current cell has DoFs at all, and if not
+    // then there is nothing else to do.
     const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
-                                             shape_hessians[0].size() : 0;
+    if (dofs_per_cell == 0)
+      return;
+
+
+    const unsigned int n_quadrature_points = shape_hessians[0].size();
     const unsigned int n_components = fe.n_components();
 
     // Assert that we can write all components into the result vectors
@@ -2505,11 +2524,6 @@ namespace internal
           AssertDimension (laplacians[i].size(), result_components);
       }
 
-    // initialize with zero
-    for (unsigned int i=0; i<laplacians.size(); ++i)
-      std::fill_n (laplacians[i].begin(), laplacians[i].size(),
-                   typename VectorType::value_type());
-
     // add up contributions of trial functions.  now check whether the shape
     // function is primitive or not. if it is, then set its only non-zero
     // component, otherwise loop over components

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.