From 4d76094bd71ec61057b161b19fef5b7657f1d1a7 Mon Sep 17 00:00:00 2001 From: heltai Date: Sun, 12 Jan 2014 10:54:07 +0000 Subject: [PATCH] Merged from trunk r32177 through r32196. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32197 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 23 +++++++- deal.II/examples/step-44/step-44.cc | 2 +- deal.II/include/deal.II/fe/fe.h | 28 ++++----- deal.II/include/deal.II/meshworker/simple.h | 38 +++++------- deal.II/source/distributed/grid_refinement.cc | 38 +++++------- deal.II/source/fe/fe.cc | 31 +++++----- deal.II/source/fe/fe_poly_tensor.cc | 12 ++-- deal.II/source/fe/fe_system.cc | 14 ++++- deal.II/source/fe/fe_values.cc | 58 ++++++++++++------- 9 files changed, 135 insertions(+), 109 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 8b705acd64..4c21c7596b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -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.
    +
  1. 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. +
    + (Wolfgang Bangerth, 2014/01/08) +
    2. + +
  2. Fixed: DataOut got confused in some situations where one uses FE_Nothing. + This is now fixed. +
    + (Minh Do-Quang, Wolfgang Bangerth, 2014/01/08) +
    3. + +
  3. 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. +
    + (Minh Do-Quang, Wolfgang Bangerth, 2014/01/08) +
    4. +
  4. 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 diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc index 5e27647810..f7b59a2f17 100644 --- a/deal.II/examples/step-44/step-44.cc +++ b/deal.II/examples/step-44/step-44.cc @@ -15,7 +15,7 @@ * * --------------------------------------------------------------------- -/* + * * Authors: Jean-Paul Pelteret, University of Cape Town, * Andrew McBride, University of Erlangen-Nuremberg, 2010 */ diff --git a/deal.II/include/deal.II/fe/fe.h b/deal.II/include/deal.II/fe/fe.h index 40b4c742c5..fda2a7fef3 100644 --- a/deal.II/include/deal.II/fe/fe.h +++ b/deal.II/include/deal.II/fe/fe.h @@ -931,6 +931,19 @@ public: std::pair 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 diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h index 224729da48..ed692c57f3 100644 --- a/deal.II/include/deal.II/meshworker/simple.h +++ b/deal.II/include/deal.II/meshworker/simple.h @@ -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)); } } diff --git a/deal.II/source/distributed/grid_refinement.cc b/deal.II/source/distributed/grid_refinement.cc index 9b89144f91..f7f185ff4e 100644 --- a/deal.II/source/distributed/grid_refinement.cc +++ b/deal.II/source/distributed/grid_refinement.cc @@ -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 number diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc index f8c0482261..43db5378c1 100644 --- a/deal.II/source/fe/fe.cc +++ b/deal.II/source/fe/fe.cc @@ -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::faces_per_cell, ExcIndexRange(face, 0, GeometryInfo::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::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::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 unsigned int FiniteElement::adjust_quad_dof_index_for_face_orientation (const unsigned int index, @@ -1248,6 +1248,11 @@ FiniteElement::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(), diff --git a/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/source/fe/fe_poly_tensor.cc index 9f6460c363..4735a210ef 100644 --- a/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/source/fe/fe_poly_tensor.cc @@ -585,9 +585,9 @@ FE_PolyTensor::fill_fe_values ( } if (flags & update_hessians && cell_similarity != CellSimilarity::translation) - this->compute_2nd (mapping, cell, - typename QProjector::DataSetDescriptor().cell(), - mapping_data, fe_data, data); + this->compute_2nd (mapping, cell, + typename QProjector::DataSetDescriptor().cell(), + mapping_data, fe_data, data); } @@ -791,7 +791,7 @@ FE_PolyTensor::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::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::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); } diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc index bf70ec76a6..8683136e66 100644 --- a/deal.II/source/fe/fe_system.cc +++ b/deal.II/source/fe/fe_system.cc @@ -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::get_interpolation_matrix ( const FiniteElement &x_source_fe, FullMatrix &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)); diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc index 066c79c4c8..aaf9fe0a0a 100644 --- a/deal.II/source/fe/fe_values.cc +++ b/deal.II/source/fe/fe_values.cc @@ -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 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()); + + // 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()); - // 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 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