From: Matthias Maier Date: Sat, 2 Sep 2017 04:42:07 +0000 (-0500) Subject: intermediate X-Git-Tag: v9.0.0-rc1~1055^2~8 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0ee0baad056216abf7501b216f20f9995d3ca9f6;p=dealii.git intermediate --- diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 62f0ae0085..b55fdcec1a 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -1972,6 +1972,7 @@ public: */ bool has_generalized_support_points () const; +//FIXME /** * Return the equivalent to get_generalized_support_points(), except * for faces. @@ -1979,6 +1980,7 @@ public: const std::vector > & get_generalized_face_support_points () const; +//FIXME /** * Return whether a finite element has defined generalized support points on * faces. If the result is true, then a call to the diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index b7658cb197..cd3a3536e9 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -566,6 +566,9 @@ namespace VectorTools * @name Interpolation and projection */ //@{ + + + /** * Compute the interpolation of @p function at the support points to the * finite element space described by the Triangulation and FiniteElement diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 9ff7b90ff5..27adfdd4fd 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -85,215 +85,282 @@ DEAL_II_NAMESPACE_OPEN namespace VectorTools { - template class DoFHandlerType> - void interpolate (const Mapping &mapping, - const DoFHandlerType &dof, - const Function &function, - VectorType &vec, - const ComponentMask &component_mask) + void interpolate( + const Mapping &mapping, + const DoFHandlerType &dof_handler, + const Function &function, + VectorType &vec, + const ComponentMask &component_mask) { - typedef typename VectorType::value_type number; - Assert (component_mask.represents_n_components(dof.get_fe(0).n_components()), - ExcMessage("The number of components in the mask has to be either " - "zero or equal to the number of components in the finite " - "element.") ); + Assert(component_mask.represents_n_components(dof_handler.get_fe().n_components()), + ExcMessage("The number of components in the mask has to be either " + "zero or equal to the number of components in the finite " + "element.")); - Assert (vec.size() == dof.n_dofs(), - ExcDimensionMismatch (vec.size(), dof.n_dofs())); - Assert (dof.get_fe(0).n_components() == function.n_components, - ExcDimensionMismatch(dof.get_fe(0).n_components(), - function.n_components)); - Assert (component_mask.n_selected_components(dof.get_fe(0).n_components()) > 0, - ComponentMask::ExcNoComponentSelected()); + Assert (vec.size() == dof_handler.n_dofs(), + ExcDimensionMismatch (vec.size(), dof_handler.n_dofs())); - const hp::FECollection &fe = dof.get_fe_collection(); - const unsigned int n_components = fe.n_components(); - const bool fe_is_system = (n_components != 1); + Assert (dof_handler.get_fe().n_components() == function.n_components, + ExcDimensionMismatch(dof_handler.get_fe().n_components(), function.n_components)); - typename DoFHandlerType::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); + Assert (component_mask.n_selected_components(dof_handler.get_fe().n_components()) > 0, + ComponentMask::ExcNoComponentSelected()); - // For FESystems many of the - // unit_support_points will appear - // multiple times, as a point may be - // unit_support_point for several of the - // components of the system. The following - // is rather complicated, but at least - // attempts to avoid evaluating the - // vectorfunction multiple times at the - // same point on a cell. // - // First check that the desired components are interpolating. + // Computing the generalized interpolant isn't quite as straightforward + // as for classical Lagrange elements. A major complication is the fact + // it generally doesn't hold true that a function evaluates to the same + // dof coefficient on different cells. This means *setting* the value + // of a (global) degree of freedom computed on one cell doesn't + // necessarily lead to the same result when computed on a neighboring + // cell (that shares the same global degree of freedom). + // + // We thus, do the following operation: + // + // On each cell: + // + // - We first determine all function values u(x_i) in generalized + // support points + // + // - We transform these function values back to the unit cell + // according to the conformity of the component (scalar, Hdiv, or + // Hcurl conforming); see [Monk, Finite Element Methods for Maxwell's + // Equations, p.77ff Section 3.9] for details. This results in + // \hat u(\hat x_i) + // + // - We convert these generalized support point values to nodal values + // + // - For every global dof we take the average 1 / n_K \sum_{K} dof_K + // where n_K is the number of cells sharing the global dof and dof_K + // is the computed value on the cell K. + // + // For every degree of freedom that is shared by k cells, we compute + // its value on all k cells and take the weighted average with respect + // to the JxW values. + // + + typedef typename VectorType::value_type number; + + const hp::FECollection fe(dof_handler.get_fe()); + + std::vector dofs_on_cell(fe.max_dofs_per_cell()); + + // Temporary storage for cell-wise interpolation operation. We store a + // variant for every fe we encounter to speed up resizing operations. + // The first vector is used for local function evaluation. The vector + // dof_values is used to store intermediate cell-wise interpolation + // results (see the detailed explanation in the for loop further down + // below). + + std::vector > > fe_function_values(fe.size()); + std::vector > fe_dof_values(fe.size()); + + // We will need two temporary global vectors that store the new values + // and weights. + VectorType interpolation; + VectorType weights; + interpolation.reinit(vec); + weights.reinit(vec); + + // We use an FEValues object to transform all generalized support + // points from the unit cell to the real cell coordinates. Thus, + // initialize a quadrature with all generalized support points and + // create an FEValues object with it. + + hp::QCollection support_quadrature; for (unsigned int fe_index=0; fe_index(points)); } - // Find the support points on a cell that are mentioned multiple times, and - // only add each once. Each multiple point gets to know the dof index of - // its representative point by the dof_to_rep_dof_table. + hp::MappingCollection mapping_collection(mapping); - // the following vector collects all unit support points p[i], - // 0<=i > > rep_unit_support_points (fe.size()); - // the following table converts a dof i - // to the rep index. - std::vector > dof_to_rep_index_table(fe.size()); + hp::FEValues fe_values( + mapping_collection, + fe, + support_quadrature, + update_quadrature_points); + + // An extra FEValues object to compute jacobians. + // Only re-initialized in case of Hcurl or Hdiv conforming elements, + // i.e. if we really need the information. + hp::FEValues fe_values_jacobians( + mapping_collection, + fe, + support_quadrature, + update_jacobians | update_inverse_jacobians); - std::vector n_rep_points (fe.size(), 0); + // + // Now loop over all locally owned, active cells. + // - for (unsigned int fe_index=0; fe_indexis_locally_owned()) + continue; + + const unsigned int fe_index = cell->active_fe_index(); + + // Do nothing if there are no local degrees of freedom. + if (fe[fe_index].dofs_per_cell == 0) + continue; + + // Get transformed, generalized support points + fe_values.reinit(cell); + const std::vector > &generalized_support_points = + fe_values.get_present_fe_values().get_quadrature_points(); + + // Get indices of the dofs on this cell + const auto n_dofs = fe[fe_index].dofs_per_cell; + dofs_on_cell.resize (n_dofs); + cell->get_dof_indices (dofs_on_cell); + + // Prepare temporary storage + auto &function_values = fe_function_values[fe_index]; + auto &dof_values = fe_dof_values[fe_index]; + + const auto n_components = fe[fe_index].n_components(); + function_values.resize(generalized_support_points.size(), + Vector(n_components)); + dof_values.resize(n_dofs); + + // Get all function values: + function.vector_value_list(generalized_support_points, + function_values); + + // FIXME: In case of an FESystem we have to apply this + // transformation according to the conformity of each base element. + + // Before we can average, we have to transform all function values + // from the real cell back to the unit cell. We query the finite + // element for the correct transformation. + switch (fe[fe_index].conforming_space) { - const unsigned int component - = fe[fe_index].system_to_component_index(i).first; - if (component_mask[component] == true) + case FiniteElementData::H1: + DEAL_II_FALLTHROUGH; + case FiniteElementData::L2: + // See Monk, Finite Element Methods for Maxwell's Equations, + // p. 77ff, formula (3.74). + // For given mapping F_K: \hat K \to K, we have to transform + // \hat p = p\circ F_K + // i.e., do nothing. + break; + + case FiniteElementData::Hcurl: + // See Monk, Finite Element Methods for Maxwell's Equations, + // p. 77ff, formula (3.76) and Corollary 3.58. + // For given mapping F_K: \hat K \to K, we have to transform + // \hat u = (dF_K)^T u\circ F_K + + fe_values_jacobians.reinit(cell); + for (unsigned int i = 0; i < function_values.size(); ++i) { - const Point dof_support_point = fe[fe_index].unit_support_point(i); - - bool representative=true; - // the following loop is looped - // the other way round to get - // the minimal effort of - // O(fe.dofs_per_cell) for multiple - // support points that are placed - // one after the other. - for (unsigned int j=rep_unit_support_points[fe_index].size(); j>0; --j) - if (dof_support_point - == rep_unit_support_points[fe_index][j-1]) - { - dof_to_rep_index_table[fe_index].push_back(j-1); - representative=false; - break; - } + const auto &jacobians = + fe_values_jacobians.get_present_fe_values().get_jacobians(); - if (representative) - { - dof_to_rep_index_table[fe_index].push_back - (rep_unit_support_points[fe_index].size()); - rep_unit_support_points[fe_index].push_back(dof_support_point); - ++n_rep_points[fe_index]; - } + // value[m] <- sum jacobian_transpose[m][n] * old_value[n]: + const auto old_value = function_values[i]; + TensorAccessors::contract<1, 2, 1, dim>( + function_values[i], jacobians[i].transpose(), old_value); } - else + break; + + case FiniteElementData::Hdiv: + // See Monk, Finite Element Methods for Maxwell's Equations, + // p. 79ff, formula (3.77) and Lemma 3.59. + // For given mapping F_K: \hat K \to K, we have to transform + // \hat w = det(dF_K) (dF_K)^{-1} w\circ F_K + + fe_values_jacobians.reinit(cell); + for (unsigned int i = 0; i < function_values.size(); ++i) { - // If correct component not to be interpolated, append invalid index to in table - dof_to_rep_index_table[fe_index].push_back(numbers::invalid_dof_index); + const auto &jacobians = + fe_values_jacobians.get_present_fe_values().get_jacobians(); + const auto &inverse_jacobians = + fe_values_jacobians.get_present_fe_values().get_inverse_jacobians(); + + // value[m] <- sum inverse_jacobians[m][n] * old_value[n]: + const auto old_value = function_values[i]; + TensorAccessors::contract<1, 2, 1, dim>( + function_values[i], inverse_jacobians[i], old_value); + + for (unsigned int j = 0; j < n_components; ++j) + function_values[i][j] *= jacobians[i].determinant(); } - } - - Assert(rep_unit_support_points[fe_index].size()==n_rep_points[fe_index], - ExcInternalError()); - Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell, - ExcInternalError()); - } + break; - std::vector dofs_on_cell (fe.max_dofs_per_cell()); + default: + Assert(false, + ExcMessage( + "The supplied finite element has an unknown conformity.")); + } /*switch*/ - // get space for the values of the - // function at the rep support points. - // - // have two versions, one for system fe - // and one for scalar ones, to take the - // more efficient one respectively - std::vector > function_values_scalar(fe.size()); - std::vector > > function_values_system(fe.size()); - - // Make a quadrature rule from support points - // to feed it into FEValues - hp::QCollection support_quadrature; - for (unsigned int fe_index=0; fe_index(rep_unit_support_points[fe_index])); + FETools::convert_generalized_support_point_values_to_dof_values( + fe[fe_index], function_values, dof_values); - // Transformed support points are computed by - // FEValues - hp::MappingCollection mapping_collection (mapping); + for (unsigned int i=0; i < n_dofs; ++i) + { + ::dealii::internal::ElementAccess::add( + typename VectorType::value_type(1.0), + dofs_on_cell[i], + weights); - hp::FEValues fe_values (mapping_collection, - fe, support_quadrature, update_quadrature_points); + const auto &nonzero_components = + fe[fe_index].get_nonzero_components(i); - for (; cell!=endc; ++cell) - if (cell->is_locally_owned()) - { - const unsigned int fe_index = cell->active_fe_index(); - if (fe[fe_index].dofs_per_cell != 0) - { - // for each cell: - // get location of finite element - // support_points - fe_values.reinit(cell); - const std::vector > &rep_support_points = - fe_values.get_present_fe_values().get_quadrature_points(); + // Figure out whether the component mask applies. We assume + // that we are allowed to set degrees of freedom if at least + // one of the components (of the dof) is selected. + // + // FIXME: Add a debug assert. + bool selected = false; + for (unsigned int i = 0; i < nonzero_components.size(); ++i) + selected = + selected || (nonzero_components[i] && component_mask[i]); - // get indices of the dofs on this cell - dofs_on_cell.resize (fe[fe_index].dofs_per_cell); - cell->get_dof_indices (dofs_on_cell); + if (selected) + { + // Add local values to the global vectors + ::dealii::internal::ElementAccess::add( + dof_values[i], dofs_on_cell[i], interpolation); + } + else + { + // If a component is ignored, simply copy all dof values + // from the vector "vec": + const auto value = + ::dealii::internal::ElementAccess::get( + vec, dofs_on_cell[i]); + ::dealii::internal::ElementAccess::add( + value, dofs_on_cell[i], interpolation); + } + } + } /* loop over dof_handler.active_cell_iterators() */ + interpolation.compress(VectorOperation::add); + weights.compress(VectorOperation::add); - if (fe_is_system) - { - // get function values at - // these points. Here: get - // all components - function_values_system[fe_index] - .resize (n_rep_points[fe_index], - Vector (fe[fe_index].n_components())); - function.vector_value_list (rep_support_points, - function_values_system[fe_index]); - // distribute the function - // values to the global - // vector - for (unsigned int i=0; i::set( - function_values_system[fe_index][rep_dof](component), - dofs_on_cell[i], vec); - } - } - } - else - { - // get first component only, - // which is the only component - // in the function anyway - function_values_scalar[fe_index].resize (n_rep_points[fe_index]); - function.value_list (rep_support_points, - function_values_scalar[fe_index], - 0); - // distribute the function - // values to the global - // vector - for (unsigned int i=0; i::set( - function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]], - dofs_on_cell[i], vec); - } - } - } + for (const auto i : interpolation.locally_owned_elements()) + { + const auto value = + ::dealii::internal::ElementAccess::get(interpolation, i); + const auto weight = + ::dealii::internal::ElementAccess::get(weights, i); + ::dealii::internal::ElementAccess::set( + value / weight, i, vec); + } vec.compress(VectorOperation::insert); } + + template class DoFHandlerType> void interpolate (const DoFHandlerType &dof, diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in index 6e1d4d615a..b8e1bf7964 100644 --- a/source/numerics/vector_tools_interpolate.inst.in +++ b/source/numerics/vector_tools_interpolate.inst.in @@ -14,6 +14,7 @@ // --------------------------------------------------------------------- + for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension