From 76082c884a58197e0ba8a6c4362448a06bf1c4c3 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Thu, 7 Sep 2017 13:10:34 -0500 Subject: [PATCH] Refactor template definitions into fe_tools_templates.h header --- include/deal.II/fe/fe_tools.h | 6 +- include/deal.II/fe/fe_tools.templates.h | 3130 +++++++++++++++++++++++ source/fe/fe_tools.cc | 3110 +--------------------- 3 files changed, 3132 insertions(+), 3114 deletions(-) create mode 100644 include/deal.II/fe/fe_tools.templates.h diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index 8e8228eda1..95bb7ff008 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -547,7 +547,6 @@ namespace FETools - /** * This method implements the * FETools::compute_projection_from_quadrature_points_matrix method for @@ -1351,7 +1350,4 @@ namespace FETools DEAL_II_NAMESPACE_CLOSE -/*---------------------------- fe_tools.h ---------------------------*/ -/* end of #ifndef dealii_fe_tools_H */ -#endif -/*---------------------------- fe_tools.h ---------------------------*/ +#endif /* dealii_fe_tools_H */ diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h new file mode 100644 index 0000000000..7bed398ebc --- /dev/null +++ b/include/deal.II/fe/fe_tools.templates.h @@ -0,0 +1,3130 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_fe_tools_templates_H +#define dealii_fe_tools_templates_H + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include + +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace FETools +{ + namespace Compositing + { + + template + FiniteElementData + multiply_dof_numbers (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product) + { + AssertDimension(fes.size(), multiplicities.size()); + + unsigned int multiplied_dofs_per_vertex = 0; + unsigned int multiplied_dofs_per_line = 0; + unsigned int multiplied_dofs_per_quad = 0; + unsigned int multiplied_dofs_per_hex = 0; + + unsigned int multiplied_n_components = 0; + + unsigned int degree = 0; // degree is the maximal degree of the components + + unsigned int n_components = 0; + // Get the number of components from the first given finite element. + for (unsigned int i=0; i0) + { + n_components = fes[i]->n_components(); + break; + } + + for (unsigned int i=0; i0) + { + multiplied_dofs_per_vertex += fes[i]->dofs_per_vertex * multiplicities[i]; + multiplied_dofs_per_line += fes[i]->dofs_per_line * multiplicities[i]; + multiplied_dofs_per_quad += fes[i]->dofs_per_quad * multiplicities[i]; + multiplied_dofs_per_hex += fes[i]->dofs_per_hex * multiplicities[i]; + + multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; + + Assert (do_tensor_product || (n_components == fes[i]->n_components()), + ExcDimensionMismatch(n_components, fes[i]->n_components())); + + degree = std::max(degree, fes[i]->tensor_degree() ); + } + + // assume conformity of the first finite element and then take away + // bits as indicated by the base elements. if all multiplicities + // happen to be zero, then it doesn't matter what we set it to. + typename FiniteElementData::Conformity total_conformity + = typename FiniteElementData::Conformity(); + { + unsigned int index = 0; + for (index=0; index0) + { + total_conformity = fes[index]->conforming_space; + break; + } + + for (; index0) + total_conformity = + typename FiniteElementData::Conformity(total_conformity + & + fes[index]->conforming_space); + } + + std::vector dpo; + dpo.push_back(multiplied_dofs_per_vertex); + dpo.push_back(multiplied_dofs_per_line); + if (dim>1) dpo.push_back(multiplied_dofs_per_quad); + if (dim>2) dpo.push_back(multiplied_dofs_per_hex); + + BlockIndices block_indices (0,0); + + for (unsigned int base=0; base < fes.size(); ++base) + for (unsigned int m = 0; m < multiplicities[base]; ++m) + block_indices.push_back(fes[base]->dofs_per_cell); + + return FiniteElementData (dpo, + (do_tensor_product ? multiplied_n_components : n_components), + degree, + total_conformity, + block_indices); + } + + + + template + FiniteElementData + multiply_dof_numbers (const FiniteElement *fe1, + const unsigned int N1, + const FiniteElement *fe2, + const unsigned int N2, + const FiniteElement *fe3, + const unsigned int N3, + const FiniteElement *fe4, + const unsigned int N4, + const FiniteElement *fe5, + const unsigned int N5) + { + std::vector*> fes; + fes.push_back(fe1); + fes.push_back(fe2); + fes.push_back(fe3); + fes.push_back(fe4); + fes.push_back(fe5); + + std::vector mult; + mult.push_back(N1); + mult.push_back(N2); + mult.push_back(N3); + mult.push_back(N4); + mult.push_back(N5); + return multiply_dof_numbers(fes, mult); + } + + + + template + std::vector + compute_restriction_is_additive_flags (const std::vector*> &fes, + const std::vector &multiplicities) + { + AssertDimension(fes.size(), multiplicities.size()); + + // first count the number of dofs and components that will emerge from the + // given FEs + unsigned int n_shape_functions = 0; + for (unsigned int i=0; i0) // check needed as fe might be NULL + n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; + + // generate the array that will hold the output + std::vector retval (n_shape_functions, false); + + // finally go through all the shape functions of the base elements, and copy + // their flags. this somehow copies the code in build_cell_table, which is + // not nice as it uses too much implicit knowledge about the layout of the + // individual bases in the composed FE, but there seems no way around... + // + // for each shape function, copy the flags from the base element to this + // one, taking into account multiplicities, and other complications + unsigned int total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; + ++vertex_number) + { + for (unsigned int base=0; basedofs_per_vertex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 2. Lines + if (GeometryInfo::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + for (unsigned int base=0; basedofs_per_line; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_line*line_number + + local_index + + fes[base]->first_line_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 3. Quads + if (GeometryInfo::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + for (unsigned int base=0; basedofs_per_quad; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_quad*quad_number + + local_index + + fes[base]->first_quad_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + // 4. Hexes + if (GeometryInfo::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + for (unsigned int base=0; basedofs_per_hex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_hex*hex_number + + local_index + + fes[base]->first_hex_index); + + Assert (index_in_base < fes[base]->dofs_per_cell, + ExcInternalError()); + retval[total_index] = fes[base]->restriction_is_additive(index_in_base); + } + } + + Assert (total_index == n_shape_functions, ExcInternalError()); + + return retval; + } + + + + /** + * Take a @p FiniteElement object + * and return an boolean vector including the @p + * restriction_is_additive_flags of the mixed element consisting of @p N + * elements of the sub-element @p fe. + */ + template + std::vector + compute_restriction_is_additive_flags (const FiniteElement *fe1, + const unsigned int N1, + const FiniteElement *fe2, + const unsigned int N2, + const FiniteElement *fe3, + const unsigned int N3, + const FiniteElement *fe4, + const unsigned int N4, + const FiniteElement *fe5, + const unsigned int N5) + { + std::vector*> fe_list; + std::vector multiplicities; + + fe_list.push_back (fe1); + multiplicities.push_back (N1); + + fe_list.push_back (fe2); + multiplicities.push_back (N2); + + fe_list.push_back (fe3); + multiplicities.push_back (N3); + + fe_list.push_back (fe4); + multiplicities.push_back (N4); + + fe_list.push_back (fe5); + multiplicities.push_back (N5); + return compute_restriction_is_additive_flags (fe_list, multiplicities); + } + + + + template + std::vector + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities, + const bool do_tensor_product) + { + AssertDimension(fes.size(), multiplicities.size()); + + // first count the number of dofs and components that will emerge from the + // given FEs + unsigned int n_shape_functions = 0; + for (unsigned int i=0; i0) //needed because fe might be NULL + n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; + + unsigned int n_components = 0; + if (do_tensor_product) + { + for (unsigned int i=0; i0) //needed because fe might be NULL + n_components += fes[i]->n_components() * multiplicities[i]; + } + else + { + for (unsigned int i=0; i0) //needed because fe might be NULL + { + n_components = fes[i]->n_components(); + break; + } + // Now check that all FEs have the same number of components: + for (unsigned int i=0; i0) //needed because fe might be NULL + Assert (n_components == fes[i]->n_components(), + ExcDimensionMismatch(n_components,fes[i]->n_components())); + } + + // generate the array that will hold the output + std::vector > + retval (n_shape_functions, std::vector (n_components, false)); + + // finally go through all the shape functions of the base elements, and copy + // their flags. this somehow copies the code in build_cell_table, which is + // not nice as it uses too much implicit knowledge about the layout of the + // individual bases in the composed FE, but there seems no way around... + // + // for each shape function, copy the non-zero flags from the base element to + // this one, taking into account multiplicities, multiple components in base + // elements, and other complications + unsigned int total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; + ++vertex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; basen_components() * do_tensor_product) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_vertex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + // 2. Lines + if (GeometryInfo::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; basen_components() * do_tensor_product) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_line; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_line*line_number + + local_index + + fes[base]->first_line_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + // 3. Quads + if (GeometryInfo::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; basen_components() * do_tensor_product) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_quad; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_quad*quad_number + + local_index + + fes[base]->first_quad_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + // 4. Hexes + if (GeometryInfo::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; basen_components() * do_tensor_product) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_hex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_hex*hex_number + + local_index + + fes[base]->first_hex_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->get_nonzero_components(index_in_base)[c]; + } + } + } + + Assert (total_index == n_shape_functions, ExcInternalError()); + + // now copy the vector > into a vector. + // this appears complicated but we do it this way since it's just + // awkward to generate ComponentMasks directly and so we need the + // recourse of the inner vector anyway. + std::vector xretval (retval.size()); + for (unsigned int i=0; i + std::vector + compute_nonzero_components (const FiniteElement *fe1, + const unsigned int N1, + const FiniteElement *fe2, + const unsigned int N2, + const FiniteElement *fe3, + const unsigned int N3, + const FiniteElement *fe4, + const unsigned int N4, + const FiniteElement *fe5, + const unsigned int N5, + const bool do_tensor_product) + { + std::vector*> fe_list; + std::vector multiplicities; + + fe_list.push_back (fe1); + multiplicities.push_back (N1); + + fe_list.push_back (fe2); + multiplicities.push_back (N2); + + fe_list.push_back (fe3); + multiplicities.push_back (N3); + + fe_list.push_back (fe4); + multiplicities.push_back (N4); + + fe_list.push_back (fe5); + multiplicities.push_back (N5); + + return compute_nonzero_components (fe_list, multiplicities, + do_tensor_product); + } + + + + template + void + build_cell_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &system_to_component_table, + std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &component_to_base_table, + const FiniteElement &fe, + const bool do_tensor_product) + { + unsigned int total_index = 0; + + if (do_tensor_product) + { + for (unsigned int base=0; base < fe.n_base_elements(); ++base) + for (unsigned int m = 0; m < fe.element_multiplicity(base); ++m) + { + for (unsigned int k=0; k + non_primitive_index (numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int); + + // First enumerate vertex indices, where we first enumerate all indices on + // the first vertex in the order of the base elements, then of the second + // vertex, etc + total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; + ++vertex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base + void + build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table, + std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table, + const FiniteElement &fe, + const bool do_tensor_product) + { + // Initialize index tables. do this in the same way as done for the cell + // tables, except that we now loop over the objects of faces + + // For non-primitive shape functions, have a special invalid index + const std::pair + non_primitive_index (numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int); + + // 1. Vertices + unsigned int total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_face; + ++vertex_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::lines_per_face > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_face; + ++line_number) + { + unsigned int comp_start = 0; + for (unsigned int base = 0; base < fe.n_base_elements(); ++base) + for (unsigned int m=0; m::quads_per_face > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_face; + ++quad_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base + FiniteElement * + FEFactory::get (const Quadrature<1> &) const + { + Assert(false, ExcNotImplemented()); + return nullptr; + } + + // Specializations for FE_Q. + template <> + FiniteElement<1, 1> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q<1>(quad); + } + + template <> + FiniteElement<2, 2> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q<2>(quad); + } + + template <> + FiniteElement<3, 3> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q<3>(quad); + } + + // Specializations for FE_Q_DG0. + template <> + FiniteElement<1, 1> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q_DG0<1>(quad); + } + + template <> + FiniteElement<2, 2> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q_DG0<2>(quad); + } + + template <> + FiniteElement<3, 3> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q_DG0<3>(quad); + } + + // Specializations for FE_Q_Bubbles. + template <> + FiniteElement<1, 1> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q_Bubbles<1>(quad); + } + + template <> + FiniteElement<2, 2> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q_Bubbles<2>(quad); + } + + template <> + FiniteElement<3, 3> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_Q_Bubbles<3>(quad); + } + + // Specializations for FE_DGQArbitraryNodes. + template <> + FiniteElement<1, 1> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<1>(quad); + } + + template <> + FiniteElement<1, 2> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<1, 2>(quad); + } + + template <> + FiniteElement<1, 3> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<1, 3>(quad); + } + + template <> + FiniteElement<2, 2> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<2>(quad); + } + + template <> + FiniteElement<2, 3> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<2, 3>(quad); + } + + template <> + FiniteElement<3, 3> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<3>(quad); + } +} + +namespace +{ + // The following three functions serve to fill the maps from element + // names to elements fe_name_map below. The first one exists because + // we have finite elements which are not implemented for nonzero + // codimension. These should be transferred to the second function + // eventually. + + template + void + fill_no_codim_fe_names (std::map > &result) + { + typedef std::shared_ptr FEFactoryPointer; + + result["FE_Q_Hierarchical"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_ABF"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Bernstein"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_BDM"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGBDM"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGNedelec"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGRaviartThomas"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_RaviartThomas"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_RaviartThomasNodal"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Nedelec"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGPNonparametric"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGP"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGPMonomial"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQ"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQArbitraryNodes"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQLegendre"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQHermite"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_FaceQ"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_FaceP"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q_DG0"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q_Bubbles"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q_iso_Q1"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Nothing"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_RannacherTurek"] + = FEFactoryPointer(new FETools::FEFactory >); + } + + + + // This function fills a map from names to finite elements for any + // dimension and codimension for those elements which support + // nonzero codimension. + template + void + fill_codim_fe_names (std::map > &result) + { + typedef std::shared_ptr FEFactoryPointer; + + result["FE_Bernstein"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGP"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQ"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Nothing"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQArbitraryNodes"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQLegendre"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_DGQHermite"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q_Bubbles"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q_DG0"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q_iso_Q1"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Q"] + = FEFactoryPointer(new FETools::FEFactory >); + result["FE_Bernstein"] + = FEFactoryPointer(new FETools::FEFactory >); + } + + // The function filling the vector fe_name_map below. It iterates + // through all legal dimension/spacedimension pairs and fills + // fe_name_map[dimension][spacedimension] with the maps generated + // by the functions above. + std::vector > > > + fill_default_map() + { + std::vector > > > + result(4); + + for (unsigned int d=0; d<4; ++d) + result[d].resize(4); + + fill_no_codim_fe_names<1> (result[1][1]); + fill_no_codim_fe_names<2> (result[2][2]); + fill_no_codim_fe_names<3> (result[3][3]); + + fill_codim_fe_names<1,2> (result[1][2]); + fill_codim_fe_names<1,3> (result[1][3]); + fill_codim_fe_names<2,3> (result[2][3]); + + return result; + } + + + // have a lock that guarantees that at most one thread is changing + // and accessing the fe_name_map variable. make this lock local to + // this file. + // + // this and the next variable are declared static (even though + // they're in an anonymous namespace) in order to make icc happy + // (which otherwise reports a multiply defined symbol when linking + // libraries for more than one space dimension together + static + Threads::Mutex fe_name_map_lock; + + // This is the map used by FETools::get_fe_by_name and + // FETools::add_fe_name. It is only accessed by functions in this + // file, so it is safe to make it a static variable here. It must be + // static so that we can link several dimensions together. + + // The organization of this storage is such that + // fe_name_map[dim][spacedim][name] points to an + // FEFactoryBase with the name given. Since + // all entries of this vector are of different type, we store + // pointers to generic objects and cast them when needed. + + // We use a shared pointer to factory objects, to ensure that they + // get deleted at the end of the program run and don't end up as + // apparent memory leaks to programs like valgrind. + + // This vector is initialized at program start time using the + // function above. because at this time there are no threads + // running, there are no thread-safety issues here. since this is + // compiled for all dimensions at once, need to create objects for + // each dimension and then separate between them further down + static + std::vector > > > + fe_name_map = fill_default_map(); +} + + + + + + +namespace +{ + + // forwarder function for + // FE::get_interpolation_matrix. we + // will want to call that function + // for arbitrary FullMatrix + // types, but it only accepts + // double arguments. since it is a + // virtual function, this can also + // not be changed. so have a + // forwarder function that calls + // that function directly if + // T==double, and otherwise uses a + // temporary + template + inline + void gim_forwarder (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + fe2.get_interpolation_matrix (fe1, interpolation_matrix); + } + + + + template + inline + void gim_forwarder (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + FullMatrix tmp (interpolation_matrix.m(), + interpolation_matrix.n()); + fe2.get_interpolation_matrix (fe1, tmp); + interpolation_matrix = tmp; + } + + + + // return how many characters + // starting at the given position + // of the string match either the + // generic string "" or the + // specialized string with "dim" + // replaced with the numeric value + // of the template argument + template + inline + unsigned int match_dimension (const std::string &name, + const unsigned int position) + { + if (position >= name.size()) + return 0; + + if ((position+5 < name.size()) + && + (name[position] == '<') + && + (name[position+1] == 'd') + && + (name[position+2] == 'i') + && + (name[position+3] == 'm') + && + (name[position+4] == '>')) + return 5; + + Assert (dim<10, ExcNotImplemented()); + const char dim_char = '0'+dim; + + if ((position+3 < name.size()) + && + (name[position] == '<') + && + (name[position+1] == dim_char) + && + (name[position+2] == '>')) + return 3; + + // some other string that doesn't + // match + return 0; + } +} + + +namespace FETools +{ + template + FEFactoryBase::~FEFactoryBase() + {} + + + + template + void compute_component_wise( + const FiniteElement &element, + std::vector &renumbering, + std::vector > &comp_start) + { + Assert(renumbering.size() == element.dofs_per_cell, + ExcDimensionMismatch(renumbering.size(), + element.dofs_per_cell)); + + comp_start.resize(element.n_base_elements()); + + unsigned int k=0; + for (unsigned int i=0; i, unsigned int> + indices = element.system_to_base_index(i); + renumbering[i] = comp_start[indices.first.first][indices.first.second] + +indices.second; + } + } + + + + template + void compute_block_renumbering ( + const FiniteElement &element, + std::vector &renumbering, + std::vector &block_data, + bool return_start_indices) + { + Assert(renumbering.size() == element.dofs_per_cell, + ExcDimensionMismatch(renumbering.size(), + element.dofs_per_cell)); + Assert(block_data.size() == element.n_blocks(), + ExcDimensionMismatch(block_data.size(), + element.n_blocks())); + + types::global_dof_index k=0; + unsigned int count=0; + for (unsigned int b=0; b start_indices(block_data.size()); + k = 0; + for (unsigned int i=0; i + indices = element.system_to_block_index(i); + renumbering[i] = start_indices[indices.first] + +indices.second; + } + } + + + + template + void get_interpolation_matrix (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(interpolation_matrix.m()==fe2.dofs_per_cell && + interpolation_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(interpolation_matrix.m(), + interpolation_matrix.n(), + fe2.dofs_per_cell, + fe1.dofs_per_cell)); + + // first try the easy way: maybe + // the FE wants to implement things + // itself: + bool fe_implements_interpolation = true; + try + { + gim_forwarder (fe1, fe2, interpolation_matrix); + } + catch (typename FiniteElement::ExcInterpolationNotImplemented &) + { + // too bad.... + fe_implements_interpolation = false; + } + if (fe_implements_interpolation == true) + return; + + // uh, so this was not the + // case. hm. then do it the hard + // way. note that this will only + // work if the element is + // primitive, so check this first + Assert (fe1.is_primitive() == true, ExcFENotPrimitive()); + Assert (fe2.is_primitive() == true, ExcFENotPrimitive()); + + // Initialize FEValues for fe1 at + // the unit support points of the + // fe2 element. + const std::vector > & + fe2_support_points = fe2.get_unit_support_points (); + + typedef FiniteElement FEL; + Assert(fe2_support_points.size()==fe2.dofs_per_cell, + typename FEL::ExcFEHasNoSupportPoints()); + + for (unsigned int i=0; i + void get_back_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(interpolation_matrix.m()==fe1.dofs_per_cell && + interpolation_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(interpolation_matrix.m(), + interpolation_matrix.n(), + fe1.dofs_per_cell, + fe1.dofs_per_cell)); + + FullMatrix first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell); + FullMatrix second_matrix(fe1.dofs_per_cell, fe2.dofs_per_cell); + + get_interpolation_matrix(fe1, fe2, first_matrix); + get_interpolation_matrix(fe2, fe1, second_matrix); + + // int_matrix=second_matrix*first_matrix + second_matrix.mmult(interpolation_matrix, first_matrix); + } + + + + template + void get_interpolation_difference_matrix (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &difference_matrix) + { + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(difference_matrix.m()==fe1.dofs_per_cell && + difference_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(difference_matrix.m(), + difference_matrix.n(), + fe1.dofs_per_cell, + fe1.dofs_per_cell)); + + FullMatrix interpolation_matrix(fe1.dofs_per_cell); + get_back_interpolation_matrix(fe1, fe2, interpolation_matrix); + + for (unsigned int i=0; i + void get_projection_matrix (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &matrix) + { + Assert (fe1.n_components() == 1, ExcNotImplemented()); + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(matrix.m()==fe2.dofs_per_cell && matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(matrix.m(), matrix.n(), + fe2.dofs_per_cell, + fe1.dofs_per_cell)); + matrix = 0; + + unsigned int n1 = fe1.dofs_per_cell; + unsigned int n2 = fe2.dofs_per_cell; + + // First, create a local mass matrix for + // the unit cell + Triangulation tr; + GridGenerator::hyper_cube(tr); + + // Choose a quadrature rule + // Gauss is exact up to degree 2n-1 + const unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree()); + Assert (degree != numbers::invalid_unsigned_int, + ExcNotImplemented()); + + QGauss quadrature(degree+1); + // Set up FEValues. + const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values; + FEValues val1 (fe1, quadrature, update_values); + val1.reinit (tr.begin_active()); + FEValues val2 (fe2, quadrature, flags); + val2.reinit (tr.begin_active()); + + // Integrate and invert mass matrix + // This happens in the target space + FullMatrix mass (n2, n2); + + for (unsigned int k=0; k b(n2); + Vector x(n2); + + for (unsigned int j=0; j + FullMatrix + compute_node_matrix(const FiniteElement &fe) + { + const unsigned int n_dofs = fe.dofs_per_cell; + + FullMatrix N (n_dofs, n_dofs); + + Assert (fe.has_generalized_support_points(), ExcNotInitialized()); + Assert (fe.n_components() == dim, ExcNotImplemented()); + + const std::vector > &points = fe.get_generalized_support_points(); + + // We need the values of the polynomials in all generalized support points. + // This function specifically works for the case where shape functions + // have 'dim' vector components, so allocate that much space + std::vector > + support_point_values (points.size(), Vector(dim)); + + // In this vector, we store the + // result of the interpolation + std::vector nodal_values(n_dofs); + + // Get the values of each shape function in turn. Remember that these + // are the 'raw' shape functions (i.e., where the element has not yet + // computed the expansion coefficients with regard to the basis + // provided by the polynomial space). + for (unsigned int i=0; i + void + compute_embedding_matrices(const FiniteElement<1,2> &, + std::vector > > &, + const bool) + { + Assert(false, ExcNotImplemented()); + } + + + template <> + void + compute_embedding_matrices(const FiniteElement<1,3> &, + std::vector > > &, + const bool) + { + Assert(false, ExcNotImplemented()); + } + + + + template <> + void + compute_embedding_matrices(const FiniteElement<2,3>&, + std::vector > >&, + const bool) + { + Assert(false, ExcNotImplemented()); + } + + */ + + namespace + { + template + void + compute_embedding_for_shape_function ( + const unsigned int i, + const FiniteElement &fe, + const FEValues &coarse, + const Householder &H, + FullMatrix &this_matrix, + const double threshold) + { + const unsigned int n = fe.dofs_per_cell; + const unsigned int nd = fe.n_components (); + const unsigned int nq = coarse.n_quadrature_points; + + Vector v_coarse(nq*nd); + Vector v_fine(n); + + // The right hand side of + // the least squares + // problem consists of the + // function values of the + // coarse grid function in + // each quadrature point. + if (fe.is_primitive ()) + { + const unsigned int + d = fe.system_to_component_index (i).first; + const double *phi_i = &coarse.shape_value (i, 0); + + for (unsigned int k = 0; k < nq; ++k) + v_coarse (k * nd + d) = phi_i[k]; + } + + else + for (unsigned int d = 0; d < nd; ++d) + for (unsigned int k = 0; k < nq; ++k) + v_coarse (k * nd + d) = coarse.shape_value_component (i, k, d); + + // solve the least squares + // problem. + const double result = H.least_squares (v_fine, v_coarse); + Assert (result <= threshold, ExcLeastSquaresError (result)); + // Avoid warnings in release mode + (void)result; + (void)threshold; + + // Copy into the result + // matrix. Since the matrix + // maps a coarse grid + // function to a fine grid + // function, the columns + // are fine grid. + for (unsigned int j = 0; j < n; ++j) + this_matrix(j, i) = v_fine(j); + } + + + + template + void + compute_embedding_matrices_for_refinement_case ( + const FiniteElement &fe, + std::vector > &matrices, + const unsigned int ref_case, + const double threshold) + { + const unsigned int n = fe.dofs_per_cell; + const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); + for (unsigned int i = 0; i < nc; ++i) + { + Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n (), n)); + Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m (), n)); + } + + // Set up meshes, one with a single + // reference cell and refine it once + Triangulation tria; + GridGenerator::hyper_cube (tria, 0, 1); + tria.begin_active()->set_refine_flag (RefinementCase(ref_case)); + tria.execute_coarsening_and_refinement (); + + const unsigned int degree = fe.degree; + QGauss q_fine (degree+1); + const unsigned int nq = q_fine.size(); + + FEValues fine (fe, q_fine, + update_quadrature_points | + update_JxW_values | + update_values); + + // We search for the polynomial on + // the small cell, being equal to + // the coarse polynomial in all + // quadrature points. + + // First build the matrix for this + // least squares problem. This + // contains the values of the fine + // cell polynomials in the fine + // cell grid points. + + // This matrix is the same for all + // children. + fine.reinit (tria.begin_active ()); + const unsigned int nd = fe.n_components (); + FullMatrix A (nq*nd, n); + + for (unsigned int j = 0; j < n; ++j) + for (unsigned int d = 0; d < nd; ++d) + for (unsigned int k = 0; k < nq; ++k) + A (k * nd + d, j) = fine.shape_value_component (j, k, d); + + Householder H (A); + unsigned int cell_number = 0; + + Threads::TaskGroup task_group; + + for (typename Triangulation::active_cell_iterator + fine_cell = tria.begin_active (); fine_cell != tria.end (); + ++fine_cell, ++cell_number) + { + fine.reinit (fine_cell); + + // evaluate on the coarse cell (which + // is the first -- inactive -- cell on + // the lowest level of the + // triangulation we have created) + const std::vector > &q_points_fine = fine.get_quadrature_points(); + std::vector > q_points_coarse(q_points_fine.size()); + for (unsigned int i=0; i q_coarse (q_points_coarse, + fine.get_JxW_values ()); + FEValues coarse (fe, q_coarse, update_values); + + coarse.reinit (tria.begin (0)); + + FullMatrix &this_matrix = matrices[cell_number]; + + // Compute this once for each + // coarse grid basis function. can + // spawn subtasks if n is + // sufficiently large so that there + // are more than about 5000 + // operations in the inner loop + // (which is basically const * n^2 + // operations). + if (n > 30) + { + for (unsigned int i = 0; i < n; ++i) + { + task_group += + Threads::new_task (&compute_embedding_for_shape_function, + i, fe, coarse, H, this_matrix, threshold); + } + task_group.join_all(); + } + else + { + for (unsigned int i = 0; i < n; ++i) + { + compute_embedding_for_shape_function + (i, fe, coarse, H, this_matrix, threshold); + } + } + + // Remove small entries from + // the matrix + for (unsigned int i = 0; i < this_matrix.m (); ++i) + for (unsigned int j = 0; j < this_matrix.n (); ++j) + if (std::fabs (this_matrix (i, j)) < 1e-12) + this_matrix (i, j) = 0.; + } + + Assert (cell_number == GeometryInfo::n_children (RefinementCase (ref_case)), + ExcInternalError ()); + } + } + + + + template + void + compute_embedding_matrices(const FiniteElement &fe, + std::vector > > &matrices, + const bool isotropic_only, + const double threshold) + { + Threads::TaskGroup task_group; + + // loop over all possible refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; + + for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) + task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case, + fe, matrices[ref_case-1], ref_case, threshold); + + task_group.join_all (); + } + + + + template + void + compute_face_embedding_matrices(const FiniteElement &fe, + FullMatrix (&matrices)[GeometryInfo::max_children_per_face], + const unsigned int face_coarse, + const unsigned int face_fine, + const double threshold) + { + Assert(face_coarse==0, ExcNotImplemented()); + Assert(face_fine==0, ExcNotImplemented()); + + const unsigned int nc = GeometryInfo::max_children_per_face; + const unsigned int n = fe.dofs_per_face; + const unsigned int nd = fe.n_components(); + const unsigned int degree = fe.degree; + + const bool normal = fe.conforms(FiniteElementData::Hdiv); + const bool tangential = fe.conforms(FiniteElementData::Hcurl); + + for (unsigned int i=0; i face_c_dofs(n); + std::vector face_f_dofs(n); + { + unsigned int face_dof=0; + for (unsigned int i=0; i::vertices_per_face; ++i) + { + const unsigned int offset_c = GeometryInfo::face_to_cell_vertices(face_coarse, i) + *fe.dofs_per_vertex; + const unsigned int offset_f = GeometryInfo::face_to_cell_vertices(face_fine, i) + *fe.dofs_per_vertex; + for (unsigned int j=0; j::lines_per_face; ++i) + { + const unsigned int offset_c = fe.first_line_index + + GeometryInfo::face_to_cell_lines(face_coarse, i-1) + *fe.dofs_per_line; + const unsigned int offset_f = fe.first_line_index + + GeometryInfo::face_to_cell_lines(face_fine, i-1) + *fe.dofs_per_line; + for (unsigned int j=0; j::quads_per_face; ++i) + { + const unsigned int offset_c = fe.first_quad_index + + face_coarse + *fe.dofs_per_quad; + const unsigned int offset_f = fe.first_quad_index + + face_fine + *fe.dofs_per_quad; + for (unsigned int j=0; j tria; + GridGenerator::hyper_cube (tria, 0, 1); + tria.refine_global(1); + MappingCartesian mapping; + + // Setup quadrature and FEValues + // for a face. We cannot use + // FEFaceValues and + // FESubfaceValues because of + // some nifty handling of + // refinement cases. Guido stops + // disliking and instead starts + // hating the anisotropic implementation + QGauss q_gauss(degree+1); + const Quadrature q_fine = QProjector::project_to_face(q_gauss, face_fine); + const unsigned int nq = q_fine.size(); + + FEValues fine (mapping, fe, q_fine, + update_quadrature_points | update_JxW_values | update_values); + + // We search for the polynomial on + // the small cell, being equal to + // the coarse polynomial in all + // quadrature points. + + // First build the matrix for this + // least squares problem. This + // contains the values of the fine + // cell polynomials in the fine + // cell grid points. + + // This matrix is the same for all + // children. + fine.reinit(tria.begin_active()); + FullMatrix A(nq*nd, n); + for (unsigned int j=0; j H(A); + + Vector v_coarse(nq*nd); + Vector v_fine(n); + + + + for (unsigned int cell_number = 0; cell_number < GeometryInfo::max_children_per_face; + ++cell_number) + { + const Quadrature q_coarse + = QProjector::project_to_subface(q_gauss, face_coarse, cell_number); + FEValues coarse (mapping, fe, q_coarse, update_values); + + typename Triangulation::active_cell_iterator fine_cell + = tria.begin(0)->child(GeometryInfo::child_cell_on_face( + tria.begin(0)->refinement_case(), face_coarse, cell_number)); + fine.reinit(fine_cell); + coarse.reinit(tria.begin(0)); + + FullMatrix &this_matrix = matrices[cell_number]; + + // Compute this once for each + // coarse grid basis function + for (unsigned int i=0; i + void + compute_projection_matrices(const FiniteElement &fe, + std::vector > > &matrices, + const bool isotropic_only) + { + const unsigned int n = fe.dofs_per_cell; + const unsigned int nd = fe.n_components(); + const unsigned int degree = fe.degree; + + // prepare FEValues, quadrature etc on + // coarse cell + QGauss q_fine(degree+1); + const unsigned int nq = q_fine.size(); + + // create mass matrix on coarse cell. + FullMatrix mass(n, n); + { + // set up a triangulation for coarse cell + Triangulation tr; + GridGenerator::hyper_cube (tr, 0, 1); + + FEValues coarse (fe, q_fine, + update_JxW_values | update_values); + + typename Triangulation::cell_iterator coarse_cell + = tr.begin(0); + coarse.reinit (coarse_cell); + + const std::vector &JxW = coarse.get_JxW_values(); + for (unsigned int i=0; i::isotropic_refinement + : RefinementCase::cut_x; + for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) + { + const unsigned int + nc = GeometryInfo::n_children(RefinementCase(ref_case)); + + for (unsigned int i=0; i tr; + GridGenerator::hyper_cube (tr, 0, 1); + tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); + tr.execute_coarsening_and_refinement(); + + FEValues fine (StaticMappingQ1::mapping, fe, q_fine, + update_quadrature_points | update_JxW_values | + update_values); + + typename Triangulation::cell_iterator coarse_cell + = tr.begin(0); + + Vector v_coarse(n); + Vector v_fine(n); + + for (unsigned int cell_number=0; cell_number &this_matrix = matrices[ref_case-1][cell_number]; + + // Compute right hand side, + // which is a fine level basis + // function tested with the + // coarse level functions. + fine.reinit(coarse_cell->child(cell_number)); + const std::vector > &q_points_fine = fine.get_quadrature_points(); + std::vector > q_points_coarse(q_points_fine.size()); + for (unsigned int q=0; q q_coarse (q_points_coarse, + fine.get_JxW_values()); + FEValues coarse (StaticMappingQ1::mapping, fe, q_coarse, update_values); + coarse.reinit(coarse_cell); + + // Build RHS + + const std::vector &JxW = fine.get_JxW_values(); + + // Outer loop over all fine + // grid shape functions phi_j + for (unsigned int j=0; j + void + add_fe_name(const std::string ¶meter_name, + const FEFactoryBase *factory) + { + // Erase everything after the + // actual class name + std::string name = parameter_name; + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + if (name_end < name.size()) + name.erase(name_end); + // first make sure that no other + // thread intercepts the + // operation of this function; + // for this, acquire the lock + // until we quit this function + Threads::Mutex::ScopedLock lock(fe_name_map_lock); + + Assert(fe_name_map[dim][spacedim].find(name) == fe_name_map[dim][spacedim].end(), + ExcMessage("Cannot change existing element in finite element name list")); + + // Insert the normalized name into + // the map + fe_name_map[dim][spacedim][name] = + std::shared_ptr (factory); + } + + + namespace internal + { + namespace + { + // TODO: this encapsulates the call to the + // dimension-dependent fe_name_map so that we + // have a unique interface. could be done + // smarter? + template + FiniteElement * + get_fe_by_name_ext (std::string &name, + const std::map > + &fe_name_map) + { + // Extract the name of the + // finite element class, which only + // contains characters, numbers and + // underscores. + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + const std::string name_part(name, 0, name_end); + name.erase(0, name_part.size()); + + // now things get a little more + // complicated: FESystem. it's + // more complicated, since we + // have to figure out what the + // base elements are. this can + // only be done recursively + if (name_part == "FESystem") + { + // next we have to get at the + // base elements. start with + // the first. wrap the whole + // block into try-catch to + // make sure we destroy the + // pointers we got from + // recursive calls if one of + // these calls should throw + // an exception + std::vector*> base_fes; + std::vector base_multiplicities; + try + { + // Now, just the [...] + // part should be left. + if (name.size() == 0 || name[0] != '[') + throw (std::string("Invalid first character in ") + name); + do + { + // Erase the + // leading '[' or '-' + name.erase(0,1); + // Now, the name of the + // first base element is + // first... Let's get it + base_fes.push_back (get_fe_by_name_ext (name, + fe_name_map)); + // next check whether + // FESystem placed a + // multiplicity after + // the element name + if (name[0] == '^') + { + // yes. Delete the '^' + // and read this + // multiplicity + name.erase(0,1); + + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second); + // add to length, + // including the '^' + base_multiplicities.push_back (tmp.first); + } + else + // no, so + // multiplicity is + // 1 + base_multiplicities.push_back (1); + + // so that's it for + // this base + // element. base + // elements are + // separated by '-', + // and the list is + // terminated by ']', + // so loop while the + // next character is + // '-' + } + while (name[0] == '-'); + + // so we got to the end + // of the '-' separated + // list. make sure that + // we actually had a ']' + // there + if (name.size() == 0 || name[0] != ']') + throw (std::string("Invalid first character in ") + name); + name.erase(0,1); + // just one more sanity check + Assert ((base_fes.size() == base_multiplicities.size()) + && + (base_fes.size() > 0), + ExcInternalError()); + + // ok, apparently + // everything went ok. so + // generate the composed + // element + FiniteElement *system_element = nullptr; + + // uses new FESystem constructor + // which is independent of + // the number of FEs in the system + system_element = new FESystem(base_fes, base_multiplicities); + + // now we don't need the + // list of base elements + // any more + for (unsigned int i=0; isecond.get(); + const FEFactoryBase *fef=dynamic_cast*>(ptr); + return fef->get(1); + } + else + { + // Make sure no other thread + // is just adding an element + Threads::Mutex::ScopedLock lock (fe_name_map_lock); + AssertThrow (fe_name_map.find(name_part) != fe_name_map.end(), + ExcInvalidFEName(name)); + + // Now, just the (degree) + // or (Quadrature<1>(degree+1)) + // part should be left. + if (name.size() == 0 || name[0] != '(') + throw (std::string("Invalid first character in ") + name); + name.erase(0,1); + if (name[0] != 'Q') + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second+1); + const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); + const FEFactoryBase *fef=dynamic_cast*>(ptr); + return fef->get(tmp.first); + } + else + { + unsigned int position = name.find('('); + const std::string quadrature_name(name, 0, position); + name.erase(0,position+1); + if (quadrature_name.compare("QGaussLobatto") == 0) + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + // delete "))" + name.erase(0, tmp.second+2); + const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); + const FEFactoryBase *fef=dynamic_cast*>(ptr); + return fef->get(QGaussLobatto<1>(tmp.first)); + } + else if (quadrature_name.compare("QGauss") == 0) + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + // delete "))" + name.erase(0, tmp.second+2); + const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); + const FEFactoryBase *fef=dynamic_cast*>(ptr); + return fef->get(QGauss<1>(tmp.first)); + } + else if (quadrature_name.compare("QIterated") == 0) + { + // find sub-quadrature + position = name.find('('); + const std::string subquadrature_name(name, 0, position); + AssertThrow(subquadrature_name.compare("QTrapez") == 0, + ExcNotImplemented("Could not detect quadrature of name " + subquadrature_name)); + // delete "QTrapez()," + name.erase(0,position+3); + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + // delete "))" + name.erase(0, tmp.second+2); + const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); + const FEFactoryBase *fef=dynamic_cast*>(ptr); + return fef->get(QIterated<1>(QTrapez<1>(),tmp.first)); + } + else + { + AssertThrow (false,ExcNotImplemented()); + } + } + } + + + // hm, if we have come thus far, we + // didn't know what to do with the + // string we got. so do as the docs + // say: raise an exception + AssertThrow (false, ExcInvalidFEName(name)); + + // make some compilers happy that + // do not realize that we can't get + // here after throwing + return nullptr; + } + + + + template + FiniteElement *get_fe_by_name (std::string &name) + { + return get_fe_by_name_ext (name, fe_name_map[dim][spacedim]); + } + } + } + + + + template + FiniteElement * + get_fe_from_name (const std::string ¶meter_name) + { + return get_fe_by_name (parameter_name); + } + + + + template + FiniteElement * + get_fe_by_name (const std::string ¶meter_name) + { + std::string name = Utilities::trim(parameter_name); + std::size_t index = 1; + // remove spaces that are not between two word (things that match the + // regular expression [A-Za-z0-9_]) characters. + while (2 < name.size() && index < name.size() - 1) + { + if (name[index] == ' ' && + (!(std::isalnum(name[index - 1]) || name[index - 1] == '_') || + !(std::isalnum(name[index + 1]) || name[index + 1] == '_'))) + { + name.erase(index, 1); + } + else + { + ++index; + } + } + + // Create a version of the name + // string where all template + // parameters are eliminated. + for (unsigned int pos1 = name.find('<'); + pos1 < name.size(); + pos1 = name.find('<')) + { + + const unsigned int pos2 = name.find('>'); + // If there is only a single + // character between those two, + // it should be 'd' or the number + // representing the dimension. + if (pos2-pos1 == 2) + { + const char dimchar = '0' + dim; + (void)dimchar; + if (name.at(pos1+1) != 'd') + Assert (name.at(pos1+1) == dimchar, + ExcInvalidFEDimension(name.at(pos1+1), dim)); + } + else + Assert(pos2-pos1 == 4, ExcInvalidFEName(name)); + + // If pos1==pos2, then we are + // probably at the end of the + // string + if (pos2 != pos1) + name.erase(pos1, pos2-pos1+1); + } + // Replace all occurrences of "^dim" + // by "^d" to be handled by the + // next loop + for (unsigned int pos = name.find("^dim"); + pos < name.size(); + pos = name.find("^dim")) + name.erase(pos+2, 2); + + // Replace all occurrences of "^d" + // by using the actual dimension + for (unsigned int pos = name.find("^d"); + pos < name.size(); + pos = name.find("^d")) + name.at(pos+1) = '0' + dim; + + try + { + FiniteElement *fe = internal::get_fe_by_name (name); + + // Make sure the auxiliary function + // ate up all characters of the name. + AssertThrow (name.size() == 0, + ExcInvalidFEName(parameter_name + + std::string(" extra characters after " + "end of name"))); + return fe; + } + catch (const std::string &errline) + { + AssertThrow(false, ExcInvalidFEName(parameter_name + + std::string(" at ") + + errline)); + return nullptr; + } + } + + + + template + void + compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + + // first build the matrices M and Q + // described in the documentation + FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); + + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); + M_inverse.mmult (X, Q); + + Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); + Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); + } + + + + template + void + compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (I_q.m() == quadrature.size(), + ExcMessage ("Wrong matrix size")); + Assert (I_q.n() == fe.dofs_per_cell, ExcMessage ("Wrong matrix size")); + + for (unsigned int q=0; q + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, + std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes) + { + + // check that the number columns of the projection_matrix + // matches the size of the vector_of_tensors_at_qp + Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), + ExcDimensionMismatch(projection_matrix.n_cols(), + vector_of_tensors_at_qp.size())); + + // check that the number rows of the projection_matrix + // matches the size of the vector_of_tensors_at_nodes + Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), + ExcDimensionMismatch(projection_matrix.n_rows(), + vector_of_tensors_at_nodes.size())); + + // number of support points (nodes) to project to + const unsigned int n_support_points = projection_matrix.n_rows(); + // number of quadrature points to project from + const unsigned int n_quad_points = projection_matrix.n_cols(); + + // component projected to the nodes + Vector component_at_node(n_support_points); + // component at the quadrature point + Vector component_at_qp(n_quad_points); + + for (unsigned int ii = 0; ii < dim; ++ii) + { + + component_at_qp = 0; + + // populate the vector of components at the qps + // from vector_of_tensors_at_qp + // vector_of_tensors_at_qp data is in form: + // columns: 0, 1, ..., dim + // rows: 0,1,...., n_quad_points + // so extract the ii'th column of vector_of_tensors_at_qp + for (unsigned int q = 0; q < n_quad_points; ++q) + { + component_at_qp(q) = vector_of_tensors_at_qp[q][ii]; + } + + // project from the qps -> nodes + // component_at_node = projection_matrix_u * component_at_qp + projection_matrix.vmult(component_at_node, component_at_qp); + + // rewrite the projection of the components + // back into the vector of tensors + for (unsigned int nn =0; nn + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, + std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes) + { + + // check that the number columns of the projection_matrix + // matches the size of the vector_of_tensors_at_qp + Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), + ExcDimensionMismatch(projection_matrix.n_cols(), + vector_of_tensors_at_qp.size())); + + // check that the number rows of the projection_matrix + // matches the size of the vector_of_tensors_at_nodes + Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), + ExcDimensionMismatch(projection_matrix.n_rows(), + vector_of_tensors_at_nodes.size())); + + // number of support points (nodes) + const unsigned int n_support_points = projection_matrix.n_rows(); + // number of quadrature points to project from + const unsigned int n_quad_points = projection_matrix.n_cols(); + + // number of unique entries in a symmetric second-order tensor + const unsigned int n_independent_components = + SymmetricTensor<2, dim >::n_independent_components; + + // component projected to the nodes + Vector component_at_node(n_support_points); + // component at the quadrature point + Vector component_at_qp(n_quad_points); + + // loop over the number of unique dimensions of the tensor + for (unsigned int ii = 0; ii < n_independent_components; ++ii) + { + + component_at_qp = 0; + + // row-column entry of tensor corresponding the unrolled index + TableIndices<2> row_column_index = SymmetricTensor< 2, dim >::unrolled_to_component_indices(ii); + const unsigned int row = row_column_index[0]; + const unsigned int column = row_column_index[1]; + + // populate the vector of components at the qps + // from vector_of_tensors_at_qp + // vector_of_tensors_at_qp is in form: + // columns: 0, 1, ..., n_independent_components + // rows: 0,1,...., n_quad_points + // so extract the ii'th column of vector_of_tensors_at_qp + for (unsigned int q = 0; q < n_quad_points; ++q) + { + component_at_qp(q) = (vector_of_tensors_at_qp[q])[row][column]; + } + + // project from the qps -> nodes + // component_at_node = projection_matrix_u * component_at_qp + projection_matrix.vmult(component_at_node, component_at_qp); + + // rewrite the projection of the components back into the vector of tensors + for (unsigned int nn =0; nn + void compute_projection_from_face_quadrature_points_matrix( + const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face, + FullMatrix &X) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree)); + + + + // build the matrices M and Q + // described in the documentation + FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + + { + // need an FEFaceValues object to evaluate shape function + // values on the specified face. + FEFaceValues fe_face_values (fe, lhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + + for (unsigned int i=0; i fe_face_values (fe, rhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); + + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); + M_inverse.mmult (X, Q); + + Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); + Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); + } + + + + template + void + hierarchic_to_lexicographic_numbering (unsigned int degree, std::vector &h2l) + { + // number of support points in each direction + const unsigned int n = degree+1; + + const unsigned int dofs_per_cell = Utilities::fixed_power(n); + + // Assert size maches degree + AssertDimension (h2l.size(), dofs_per_cell); + + // polynomial degree + const unsigned int dofs_per_line = degree - 1; + + // the following lines of code are somewhat odd, due to the way the + // hierarchic numbering is organized. if someone would really want to + // understand these lines, you better draw some pictures where you + // indicate the indices and orders of vertices, lines, etc, along with the + // numbers of the degrees of freedom in hierarchical and lexicographical + // order + switch (dim) + { + case 1: + { + h2l[0] = 0; + h2l[1] = dofs_per_cell-1; + for (unsigned int i=2; i + void + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, + std::vector &h2l) + { + Assert (h2l.size() == fe.dofs_per_cell, + ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); + hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); + } + + + + template + std::vector + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) + { + Assert (fe.n_components() == 1, ExcInvalidFE()); + std::vector h2l(fe.dofs_per_cell); + hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); + return (h2l); + } + + + + template + void + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, + std::vector &l2h) + { + l2h = lexicographic_to_hierarchic_numbering (fe); + } + + + + template + std::vector + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) + { + return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); + } + +} // end of namespace FETools + + +DEAL_II_NAMESPACE_CLOSE + +#endif /* dealii_fe_tools_templates_H */ diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index f3a4245f25..7f170e0320 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -14,3119 +14,11 @@ // --------------------------------------------------------------------- -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include - -#include -#include -#include - +#include DEAL_II_NAMESPACE_OPEN -namespace FETools -{ - namespace Compositing - { - - template - FiniteElementData - multiply_dof_numbers (const std::vector*> &fes, - const std::vector &multiplicities, - const bool do_tensor_product) - { - AssertDimension(fes.size(), multiplicities.size()); - - unsigned int multiplied_dofs_per_vertex = 0; - unsigned int multiplied_dofs_per_line = 0; - unsigned int multiplied_dofs_per_quad = 0; - unsigned int multiplied_dofs_per_hex = 0; - - unsigned int multiplied_n_components = 0; - - unsigned int degree = 0; // degree is the maximal degree of the components - - unsigned int n_components = 0; - // Get the number of components from the first given finite element. - for (unsigned int i=0; i0) - { - n_components = fes[i]->n_components(); - break; - } - - for (unsigned int i=0; i0) - { - multiplied_dofs_per_vertex += fes[i]->dofs_per_vertex * multiplicities[i]; - multiplied_dofs_per_line += fes[i]->dofs_per_line * multiplicities[i]; - multiplied_dofs_per_quad += fes[i]->dofs_per_quad * multiplicities[i]; - multiplied_dofs_per_hex += fes[i]->dofs_per_hex * multiplicities[i]; - - multiplied_n_components+=fes[i]->n_components() * multiplicities[i]; - - Assert (do_tensor_product || (n_components == fes[i]->n_components()), - ExcDimensionMismatch(n_components, fes[i]->n_components())); - - degree = std::max(degree, fes[i]->tensor_degree() ); - } - - // assume conformity of the first finite element and then take away - // bits as indicated by the base elements. if all multiplicities - // happen to be zero, then it doesn't matter what we set it to. - typename FiniteElementData::Conformity total_conformity - = typename FiniteElementData::Conformity(); - { - unsigned int index = 0; - for (index=0; index0) - { - total_conformity = fes[index]->conforming_space; - break; - } - - for (; index0) - total_conformity = - typename FiniteElementData::Conformity(total_conformity - & - fes[index]->conforming_space); - } - - std::vector dpo; - dpo.push_back(multiplied_dofs_per_vertex); - dpo.push_back(multiplied_dofs_per_line); - if (dim>1) dpo.push_back(multiplied_dofs_per_quad); - if (dim>2) dpo.push_back(multiplied_dofs_per_hex); - - BlockIndices block_indices (0,0); - - for (unsigned int base=0; base < fes.size(); ++base) - for (unsigned int m = 0; m < multiplicities[base]; ++m) - block_indices.push_back(fes[base]->dofs_per_cell); - - return FiniteElementData (dpo, - (do_tensor_product ? multiplied_n_components : n_components), - degree, - total_conformity, - block_indices); - } - - - - template - FiniteElementData - multiply_dof_numbers (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2, - const unsigned int N2, - const FiniteElement *fe3, - const unsigned int N3, - const FiniteElement *fe4, - const unsigned int N4, - const FiniteElement *fe5, - const unsigned int N5) - { - std::vector*> fes; - fes.push_back(fe1); - fes.push_back(fe2); - fes.push_back(fe3); - fes.push_back(fe4); - fes.push_back(fe5); - - std::vector mult; - mult.push_back(N1); - mult.push_back(N2); - mult.push_back(N3); - mult.push_back(N4); - mult.push_back(N5); - return multiply_dof_numbers(fes, mult); - } - - - - template - std::vector - compute_restriction_is_additive_flags (const std::vector*> &fes, - const std::vector &multiplicities) - { - AssertDimension(fes.size(), multiplicities.size()); - - // first count the number of dofs and components that will emerge from the - // given FEs - unsigned int n_shape_functions = 0; - for (unsigned int i=0; i0) // check needed as fe might be NULL - n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; - - // generate the array that will hold the output - std::vector retval (n_shape_functions, false); - - // finally go through all the shape functions of the base elements, and copy - // their flags. this somehow copies the code in build_cell_table, which is - // not nice as it uses too much implicit knowledge about the layout of the - // individual bases in the composed FE, but there seems no way around... - // - // for each shape function, copy the flags from the base element to this - // one, taking into account multiplicities, and other complications - unsigned int total_index = 0; - for (unsigned int vertex_number=0; - vertex_number::vertices_per_cell; - ++vertex_number) - { - for (unsigned int base=0; basedofs_per_vertex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_vertex*vertex_number + - local_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - // 2. Lines - if (GeometryInfo::lines_per_cell > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_cell; - ++line_number) - { - for (unsigned int base=0; basedofs_per_line; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_line*line_number + - local_index + - fes[base]->first_line_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - // 3. Quads - if (GeometryInfo::quads_per_cell > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_cell; - ++quad_number) - { - for (unsigned int base=0; basedofs_per_quad; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_quad*quad_number + - local_index + - fes[base]->first_quad_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - // 4. Hexes - if (GeometryInfo::hexes_per_cell > 0) - for (unsigned int hex_number= 0; - hex_number != GeometryInfo::hexes_per_cell; - ++hex_number) - { - for (unsigned int base=0; basedofs_per_hex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_hex*hex_number + - local_index + - fes[base]->first_hex_index); - - Assert (index_in_base < fes[base]->dofs_per_cell, - ExcInternalError()); - retval[total_index] = fes[base]->restriction_is_additive(index_in_base); - } - } - - Assert (total_index == n_shape_functions, ExcInternalError()); - - return retval; - } - - - - /** - * Take a @p FiniteElement object - * and return an boolean vector including the @p - * restriction_is_additive_flags of the mixed element consisting of @p N - * elements of the sub-element @p fe. - */ - template - std::vector - compute_restriction_is_additive_flags (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2, - const unsigned int N2, - const FiniteElement *fe3, - const unsigned int N3, - const FiniteElement *fe4, - const unsigned int N4, - const FiniteElement *fe5, - const unsigned int N5) - { - std::vector*> fe_list; - std::vector multiplicities; - - fe_list.push_back (fe1); - multiplicities.push_back (N1); - - fe_list.push_back (fe2); - multiplicities.push_back (N2); - - fe_list.push_back (fe3); - multiplicities.push_back (N3); - - fe_list.push_back (fe4); - multiplicities.push_back (N4); - - fe_list.push_back (fe5); - multiplicities.push_back (N5); - return compute_restriction_is_additive_flags (fe_list, multiplicities); - } - - - - template - std::vector - compute_nonzero_components (const std::vector*> &fes, - const std::vector &multiplicities, - const bool do_tensor_product) - { - AssertDimension(fes.size(), multiplicities.size()); - - // first count the number of dofs and components that will emerge from the - // given FEs - unsigned int n_shape_functions = 0; - for (unsigned int i=0; i0) //needed because fe might be NULL - n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; - - unsigned int n_components = 0; - if (do_tensor_product) - { - for (unsigned int i=0; i0) //needed because fe might be NULL - n_components += fes[i]->n_components() * multiplicities[i]; - } - else - { - for (unsigned int i=0; i0) //needed because fe might be NULL - { - n_components = fes[i]->n_components(); - break; - } - // Now check that all FEs have the same number of components: - for (unsigned int i=0; i0) //needed because fe might be NULL - Assert (n_components == fes[i]->n_components(), - ExcDimensionMismatch(n_components,fes[i]->n_components())); - } - - // generate the array that will hold the output - std::vector > - retval (n_shape_functions, std::vector (n_components, false)); - - // finally go through all the shape functions of the base elements, and copy - // their flags. this somehow copies the code in build_cell_table, which is - // not nice as it uses too much implicit knowledge about the layout of the - // individual bases in the composed FE, but there seems no way around... - // - // for each shape function, copy the non-zero flags from the base element to - // this one, taking into account multiplicities, multiple components in base - // elements, and other complications - unsigned int total_index = 0; - for (unsigned int vertex_number=0; - vertex_number::vertices_per_cell; - ++vertex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_components() * do_tensor_product) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_vertex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_vertex*vertex_number + - local_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->get_nonzero_components(index_in_base)[c]; - } - } - } - - // 2. Lines - if (GeometryInfo::lines_per_cell > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_cell; - ++line_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_components() * do_tensor_product) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_line; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_line*line_number + - local_index + - fes[base]->first_line_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->get_nonzero_components(index_in_base)[c]; - } - } - } - - // 3. Quads - if (GeometryInfo::quads_per_cell > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_cell; - ++quad_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_components() * do_tensor_product) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_quad; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_quad*quad_number + - local_index + - fes[base]->first_quad_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->get_nonzero_components(index_in_base)[c]; - } - } - } - - // 4. Hexes - if (GeometryInfo::hexes_per_cell > 0) - for (unsigned int hex_number= 0; - hex_number != GeometryInfo::hexes_per_cell; - ++hex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; basen_components() * do_tensor_product) - for (unsigned int local_index = 0; - local_index < fes[base]->dofs_per_hex; - ++local_index, ++total_index) - { - const unsigned int index_in_base - = (fes[base]->dofs_per_hex*hex_number + - local_index + - fes[base]->first_hex_index); - - Assert (comp_start+fes[base]->n_components() <= - retval[total_index].size(), - ExcInternalError()); - for (unsigned int c=0; cn_components(); ++c) - { - Assert (c < fes[base]->get_nonzero_components(index_in_base).size(), - ExcInternalError()); - retval[total_index][comp_start+c] - = fes[base]->get_nonzero_components(index_in_base)[c]; - } - } - } - - Assert (total_index == n_shape_functions, ExcInternalError()); - - // now copy the vector > into a vector. - // this appears complicated but we do it this way since it's just - // awkward to generate ComponentMasks directly and so we need the - // recourse of the inner vector anyway. - std::vector xretval (retval.size()); - for (unsigned int i=0; i - std::vector - compute_nonzero_components (const FiniteElement *fe1, - const unsigned int N1, - const FiniteElement *fe2, - const unsigned int N2, - const FiniteElement *fe3, - const unsigned int N3, - const FiniteElement *fe4, - const unsigned int N4, - const FiniteElement *fe5, - const unsigned int N5, - const bool do_tensor_product) - { - std::vector*> fe_list; - std::vector multiplicities; - - fe_list.push_back (fe1); - multiplicities.push_back (N1); - - fe_list.push_back (fe2); - multiplicities.push_back (N2); - - fe_list.push_back (fe3); - multiplicities.push_back (N3); - - fe_list.push_back (fe4); - multiplicities.push_back (N4); - - fe_list.push_back (fe5); - multiplicities.push_back (N5); - - return compute_nonzero_components (fe_list, multiplicities, - do_tensor_product); - } - - - - template - void - build_cell_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &system_to_base_table, - std::vector< std::pair< unsigned int, unsigned int > > &system_to_component_table, - std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &component_to_base_table, - const FiniteElement &fe, - const bool do_tensor_product) - { - unsigned int total_index = 0; - - if (do_tensor_product) - { - for (unsigned int base=0; base < fe.n_base_elements(); ++base) - for (unsigned int m = 0; m < fe.element_multiplicity(base); ++m) - { - for (unsigned int k=0; k - non_primitive_index (numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int); - - // First enumerate vertex indices, where we first enumerate all indices on - // the first vertex in the order of the base elements, then of the second - // vertex, etc - total_index = 0; - for (unsigned int vertex_number=0; - vertex_number::vertices_per_cell; - ++vertex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base::lines_per_cell > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_cell; - ++line_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base::quads_per_cell > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_cell; - ++quad_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base::hexes_per_cell > 0) - for (unsigned int hex_number= 0; - hex_number != GeometryInfo::hexes_per_cell; - ++hex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base - void - build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &face_system_to_base_table, - std::vector< std::pair< unsigned int, unsigned int > > &face_system_to_component_table, - const FiniteElement &fe, - const bool do_tensor_product) - { - // Initialize index tables. do this in the same way as done for the cell - // tables, except that we now loop over the objects of faces - - // For non-primitive shape functions, have a special invalid index - const std::pair - non_primitive_index (numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int); - - // 1. Vertices - unsigned int total_index = 0; - for (unsigned int vertex_number=0; - vertex_number::vertices_per_face; - ++vertex_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base::lines_per_face > 0) - for (unsigned int line_number= 0; - line_number != GeometryInfo::lines_per_face; - ++line_number) - { - unsigned int comp_start = 0; - for (unsigned int base = 0; base < fe.n_base_elements(); ++base) - for (unsigned int m=0; m::quads_per_face > 0) - for (unsigned int quad_number= 0; - quad_number != GeometryInfo::quads_per_face; - ++quad_number) - { - unsigned int comp_start = 0; - for (unsigned int base=0; base - FiniteElement * - FEFactory::get (const Quadrature<1> &) const - { - Assert(false, ExcNotImplemented()); - return nullptr; - } - - // Specializations for FE_Q. - template <> - FiniteElement<1, 1> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q<1>(quad); - } - - template <> - FiniteElement<2, 2> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q<2>(quad); - } - - template <> - FiniteElement<3, 3> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q<3>(quad); - } - - // Specializations for FE_Q_DG0. - template <> - FiniteElement<1, 1> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q_DG0<1>(quad); - } - - template <> - FiniteElement<2, 2> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q_DG0<2>(quad); - } - - template <> - FiniteElement<3, 3> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q_DG0<3>(quad); - } - - // Specializations for FE_Q_Bubbles. - template <> - FiniteElement<1, 1> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q_Bubbles<1>(quad); - } - - template <> - FiniteElement<2, 2> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q_Bubbles<2>(quad); - } - - template <> - FiniteElement<3, 3> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_Q_Bubbles<3>(quad); - } - - // Specializations for FE_DGQArbitraryNodes. - template <> - FiniteElement<1, 1> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_DGQArbitraryNodes<1>(quad); - } - - template <> - FiniteElement<1, 2> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_DGQArbitraryNodes<1, 2>(quad); - } - - template <> - FiniteElement<1, 3> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_DGQArbitraryNodes<1, 3>(quad); - } - - template <> - FiniteElement<2, 2> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_DGQArbitraryNodes<2>(quad); - } - - template <> - FiniteElement<2, 3> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_DGQArbitraryNodes<2, 3>(quad); - } - - template <> - FiniteElement<3, 3> * - FEFactory >::get (const Quadrature<1> &quad) const - { - return new FE_DGQArbitraryNodes<3>(quad); - } -} - -namespace -{ - // The following three functions serve to fill the maps from element - // names to elements fe_name_map below. The first one exists because - // we have finite elements which are not implemented for nonzero - // codimension. These should be transferred to the second function - // eventually. - - template - void - fill_no_codim_fe_names (std::map > &result) - { - typedef std::shared_ptr FEFactoryPointer; - - result["FE_Q_Hierarchical"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_ABF"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Bernstein"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_BDM"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGBDM"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGNedelec"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGRaviartThomas"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_RaviartThomas"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_RaviartThomasNodal"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Nedelec"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGPNonparametric"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGP"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGPMonomial"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQ"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQArbitraryNodes"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQLegendre"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQHermite"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_FaceQ"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_FaceP"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q_DG0"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q_Bubbles"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q_iso_Q1"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Nothing"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_RannacherTurek"] - = FEFactoryPointer(new FETools::FEFactory >); - } - - - - // This function fills a map from names to finite elements for any - // dimension and codimension for those elements which support - // nonzero codimension. - template - void - fill_codim_fe_names (std::map > &result) - { - typedef std::shared_ptr FEFactoryPointer; - - result["FE_Bernstein"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGP"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQ"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Nothing"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQArbitraryNodes"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQLegendre"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_DGQHermite"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q_Bubbles"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q_DG0"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q_iso_Q1"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Q"] - = FEFactoryPointer(new FETools::FEFactory >); - result["FE_Bernstein"] - = FEFactoryPointer(new FETools::FEFactory >); - } - - // The function filling the vector fe_name_map below. It iterates - // through all legal dimension/spacedimension pairs and fills - // fe_name_map[dimension][spacedimension] with the maps generated - // by the functions above. - std::vector > > > - fill_default_map() - { - std::vector > > > - result(4); - - for (unsigned int d=0; d<4; ++d) - result[d].resize(4); - - fill_no_codim_fe_names<1> (result[1][1]); - fill_no_codim_fe_names<2> (result[2][2]); - fill_no_codim_fe_names<3> (result[3][3]); - - fill_codim_fe_names<1,2> (result[1][2]); - fill_codim_fe_names<1,3> (result[1][3]); - fill_codim_fe_names<2,3> (result[2][3]); - - return result; - } - - - // have a lock that guarantees that at most one thread is changing - // and accessing the fe_name_map variable. make this lock local to - // this file. - // - // this and the next variable are declared static (even though - // they're in an anonymous namespace) in order to make icc happy - // (which otherwise reports a multiply defined symbol when linking - // libraries for more than one space dimension together - static - Threads::Mutex fe_name_map_lock; - - // This is the map used by FETools::get_fe_by_name and - // FETools::add_fe_name. It is only accessed by functions in this - // file, so it is safe to make it a static variable here. It must be - // static so that we can link several dimensions together. - - // The organization of this storage is such that - // fe_name_map[dim][spacedim][name] points to an - // FEFactoryBase with the name given. Since - // all entries of this vector are of different type, we store - // pointers to generic objects and cast them when needed. - - // We use a shared pointer to factory objects, to ensure that they - // get deleted at the end of the program run and don't end up as - // apparent memory leaks to programs like valgrind. - - // This vector is initialized at program start time using the - // function above. because at this time there are no threads - // running, there are no thread-safety issues here. since this is - // compiled for all dimensions at once, need to create objects for - // each dimension and then separate between them further down - static - std::vector > > > - fe_name_map = fill_default_map(); -} - - - - - - -namespace -{ - - // forwarder function for - // FE::get_interpolation_matrix. we - // will want to call that function - // for arbitrary FullMatrix - // types, but it only accepts - // double arguments. since it is a - // virtual function, this can also - // not be changed. so have a - // forwarder function that calls - // that function directly if - // T==double, and otherwise uses a - // temporary - template - inline - void gim_forwarder (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix) - { - fe2.get_interpolation_matrix (fe1, interpolation_matrix); - } - - - - template - inline - void gim_forwarder (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix) - { - FullMatrix tmp (interpolation_matrix.m(), - interpolation_matrix.n()); - fe2.get_interpolation_matrix (fe1, tmp); - interpolation_matrix = tmp; - } - - - - // return how many characters - // starting at the given position - // of the string match either the - // generic string "" or the - // specialized string with "dim" - // replaced with the numeric value - // of the template argument - template - inline - unsigned int match_dimension (const std::string &name, - const unsigned int position) - { - if (position >= name.size()) - return 0; - - if ((position+5 < name.size()) - && - (name[position] == '<') - && - (name[position+1] == 'd') - && - (name[position+2] == 'i') - && - (name[position+3] == 'm') - && - (name[position+4] == '>')) - return 5; - - Assert (dim<10, ExcNotImplemented()); - const char dim_char = '0'+dim; - - if ((position+3 < name.size()) - && - (name[position] == '<') - && - (name[position+1] == dim_char) - && - (name[position+2] == '>')) - return 3; - - // some other string that doesn't - // match - return 0; - } -} - - -namespace FETools -{ - template - FEFactoryBase::~FEFactoryBase() - {} - - - - template - void compute_component_wise( - const FiniteElement &element, - std::vector &renumbering, - std::vector > &comp_start) - { - Assert(renumbering.size() == element.dofs_per_cell, - ExcDimensionMismatch(renumbering.size(), - element.dofs_per_cell)); - - comp_start.resize(element.n_base_elements()); - - unsigned int k=0; - for (unsigned int i=0; i, unsigned int> - indices = element.system_to_base_index(i); - renumbering[i] = comp_start[indices.first.first][indices.first.second] - +indices.second; - } - } - - - - template - void compute_block_renumbering ( - const FiniteElement &element, - std::vector &renumbering, - std::vector &block_data, - bool return_start_indices) - { - Assert(renumbering.size() == element.dofs_per_cell, - ExcDimensionMismatch(renumbering.size(), - element.dofs_per_cell)); - Assert(block_data.size() == element.n_blocks(), - ExcDimensionMismatch(block_data.size(), - element.n_blocks())); - - types::global_dof_index k=0; - unsigned int count=0; - for (unsigned int b=0; b start_indices(block_data.size()); - k = 0; - for (unsigned int i=0; i - indices = element.system_to_block_index(i); - renumbering[i] = start_indices[indices.first] - +indices.second; - } - } - - - - template - void get_interpolation_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix) - { - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(interpolation_matrix.m()==fe2.dofs_per_cell && - interpolation_matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(interpolation_matrix.m(), - interpolation_matrix.n(), - fe2.dofs_per_cell, - fe1.dofs_per_cell)); - - // first try the easy way: maybe - // the FE wants to implement things - // itself: - bool fe_implements_interpolation = true; - try - { - gim_forwarder (fe1, fe2, interpolation_matrix); - } - catch (typename FiniteElement::ExcInterpolationNotImplemented &) - { - // too bad.... - fe_implements_interpolation = false; - } - if (fe_implements_interpolation == true) - return; - - // uh, so this was not the - // case. hm. then do it the hard - // way. note that this will only - // work if the element is - // primitive, so check this first - Assert (fe1.is_primitive() == true, ExcFENotPrimitive()); - Assert (fe2.is_primitive() == true, ExcFENotPrimitive()); - - // Initialize FEValues for fe1 at - // the unit support points of the - // fe2 element. - const std::vector > & - fe2_support_points = fe2.get_unit_support_points (); - - typedef FiniteElement FEL; - Assert(fe2_support_points.size()==fe2.dofs_per_cell, - typename FEL::ExcFEHasNoSupportPoints()); - - for (unsigned int i=0; i - void get_back_interpolation_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix) - { - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(interpolation_matrix.m()==fe1.dofs_per_cell && - interpolation_matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(interpolation_matrix.m(), - interpolation_matrix.n(), - fe1.dofs_per_cell, - fe1.dofs_per_cell)); - - FullMatrix first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell); - FullMatrix second_matrix(fe1.dofs_per_cell, fe2.dofs_per_cell); - - get_interpolation_matrix(fe1, fe2, first_matrix); - get_interpolation_matrix(fe2, fe1, second_matrix); - - // int_matrix=second_matrix*first_matrix - second_matrix.mmult(interpolation_matrix, first_matrix); - } - - - - template - void get_interpolation_difference_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &difference_matrix) - { - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(difference_matrix.m()==fe1.dofs_per_cell && - difference_matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(difference_matrix.m(), - difference_matrix.n(), - fe1.dofs_per_cell, - fe1.dofs_per_cell)); - - FullMatrix interpolation_matrix(fe1.dofs_per_cell); - get_back_interpolation_matrix(fe1, fe2, interpolation_matrix); - - for (unsigned int i=0; i - void get_projection_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &matrix) - { - Assert (fe1.n_components() == 1, ExcNotImplemented()); - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(matrix.m()==fe2.dofs_per_cell && matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(matrix.m(), matrix.n(), - fe2.dofs_per_cell, - fe1.dofs_per_cell)); - matrix = 0; - - unsigned int n1 = fe1.dofs_per_cell; - unsigned int n2 = fe2.dofs_per_cell; - - // First, create a local mass matrix for - // the unit cell - Triangulation tr; - GridGenerator::hyper_cube(tr); - - // Choose a quadrature rule - // Gauss is exact up to degree 2n-1 - const unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree()); - Assert (degree != numbers::invalid_unsigned_int, - ExcNotImplemented()); - - QGauss quadrature(degree+1); - // Set up FEValues. - const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values; - FEValues val1 (fe1, quadrature, update_values); - val1.reinit (tr.begin_active()); - FEValues val2 (fe2, quadrature, flags); - val2.reinit (tr.begin_active()); - - // Integrate and invert mass matrix - // This happens in the target space - FullMatrix mass (n2, n2); - - for (unsigned int k=0; k b(n2); - Vector x(n2); - - for (unsigned int j=0; j - FullMatrix - compute_node_matrix(const FiniteElement &fe) - { - const unsigned int n_dofs = fe.dofs_per_cell; - - FullMatrix N (n_dofs, n_dofs); - - Assert (fe.has_generalized_support_points(), ExcNotInitialized()); - Assert (fe.n_components() == dim, ExcNotImplemented()); - - const std::vector > &points = fe.get_generalized_support_points(); - - // We need the values of the polynomials in all generalized support points. - // This function specifically works for the case where shape functions - // have 'dim' vector components, so allocate that much space - std::vector > - support_point_values (points.size(), Vector(dim)); - - // In this vector, we store the - // result of the interpolation - std::vector nodal_values(n_dofs); - - // Get the values of each shape function in turn. Remember that these - // are the 'raw' shape functions (i.e., where the element has not yet - // computed the expansion coefficients with regard to the basis - // provided by the polynomial space). - for (unsigned int i=0; i - void - compute_embedding_matrices(const FiniteElement<1,2> &, - std::vector > > &, - const bool) - { - Assert(false, ExcNotImplemented()); - } - - - template <> - void - compute_embedding_matrices(const FiniteElement<1,3> &, - std::vector > > &, - const bool) - { - Assert(false, ExcNotImplemented()); - } - - - - template <> - void - compute_embedding_matrices(const FiniteElement<2,3>&, - std::vector > >&, - const bool) - { - Assert(false, ExcNotImplemented()); - } - - */ - - namespace - { - template - void - compute_embedding_for_shape_function ( - const unsigned int i, - const FiniteElement &fe, - const FEValues &coarse, - const Householder &H, - FullMatrix &this_matrix, - const double threshold) - { - const unsigned int n = fe.dofs_per_cell; - const unsigned int nd = fe.n_components (); - const unsigned int nq = coarse.n_quadrature_points; - - Vector v_coarse(nq*nd); - Vector v_fine(n); - - // The right hand side of - // the least squares - // problem consists of the - // function values of the - // coarse grid function in - // each quadrature point. - if (fe.is_primitive ()) - { - const unsigned int - d = fe.system_to_component_index (i).first; - const double *phi_i = &coarse.shape_value (i, 0); - - for (unsigned int k = 0; k < nq; ++k) - v_coarse (k * nd + d) = phi_i[k]; - } - - else - for (unsigned int d = 0; d < nd; ++d) - for (unsigned int k = 0; k < nq; ++k) - v_coarse (k * nd + d) = coarse.shape_value_component (i, k, d); - - // solve the least squares - // problem. - const double result = H.least_squares (v_fine, v_coarse); - Assert (result <= threshold, ExcLeastSquaresError (result)); - // Avoid warnings in release mode - (void)result; - (void)threshold; - - // Copy into the result - // matrix. Since the matrix - // maps a coarse grid - // function to a fine grid - // function, the columns - // are fine grid. - for (unsigned int j = 0; j < n; ++j) - this_matrix(j, i) = v_fine(j); - } - - - - template - void - compute_embedding_matrices_for_refinement_case ( - const FiniteElement &fe, - std::vector > &matrices, - const unsigned int ref_case, - const double threshold) - { - const unsigned int n = fe.dofs_per_cell; - const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); - for (unsigned int i = 0; i < nc; ++i) - { - Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n (), n)); - Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m (), n)); - } - - // Set up meshes, one with a single - // reference cell and refine it once - Triangulation tria; - GridGenerator::hyper_cube (tria, 0, 1); - tria.begin_active()->set_refine_flag (RefinementCase(ref_case)); - tria.execute_coarsening_and_refinement (); - - const unsigned int degree = fe.degree; - QGauss q_fine (degree+1); - const unsigned int nq = q_fine.size(); - - FEValues fine (fe, q_fine, - update_quadrature_points | - update_JxW_values | - update_values); - - // We search for the polynomial on - // the small cell, being equal to - // the coarse polynomial in all - // quadrature points. - - // First build the matrix for this - // least squares problem. This - // contains the values of the fine - // cell polynomials in the fine - // cell grid points. - - // This matrix is the same for all - // children. - fine.reinit (tria.begin_active ()); - const unsigned int nd = fe.n_components (); - FullMatrix A (nq*nd, n); - - for (unsigned int j = 0; j < n; ++j) - for (unsigned int d = 0; d < nd; ++d) - for (unsigned int k = 0; k < nq; ++k) - A (k * nd + d, j) = fine.shape_value_component (j, k, d); - - Householder H (A); - unsigned int cell_number = 0; - - Threads::TaskGroup task_group; - - for (typename Triangulation::active_cell_iterator - fine_cell = tria.begin_active (); fine_cell != tria.end (); - ++fine_cell, ++cell_number) - { - fine.reinit (fine_cell); - - // evaluate on the coarse cell (which - // is the first -- inactive -- cell on - // the lowest level of the - // triangulation we have created) - const std::vector > &q_points_fine = fine.get_quadrature_points(); - std::vector > q_points_coarse(q_points_fine.size()); - for (unsigned int i=0; i q_coarse (q_points_coarse, - fine.get_JxW_values ()); - FEValues coarse (fe, q_coarse, update_values); - - coarse.reinit (tria.begin (0)); - - FullMatrix &this_matrix = matrices[cell_number]; - - // Compute this once for each - // coarse grid basis function. can - // spawn subtasks if n is - // sufficiently large so that there - // are more than about 5000 - // operations in the inner loop - // (which is basically const * n^2 - // operations). - if (n > 30) - { - for (unsigned int i = 0; i < n; ++i) - { - task_group += - Threads::new_task (&compute_embedding_for_shape_function, - i, fe, coarse, H, this_matrix, threshold); - } - task_group.join_all(); - } - else - { - for (unsigned int i = 0; i < n; ++i) - { - compute_embedding_for_shape_function - (i, fe, coarse, H, this_matrix, threshold); - } - } - - // Remove small entries from - // the matrix - for (unsigned int i = 0; i < this_matrix.m (); ++i) - for (unsigned int j = 0; j < this_matrix.n (); ++j) - if (std::fabs (this_matrix (i, j)) < 1e-12) - this_matrix (i, j) = 0.; - } - - Assert (cell_number == GeometryInfo::n_children (RefinementCase (ref_case)), - ExcInternalError ()); - } - } - - - - template - void - compute_embedding_matrices(const FiniteElement &fe, - std::vector > > &matrices, - const bool isotropic_only, - const double threshold) - { - Threads::TaskGroup task_group; - - // loop over all possible refinement cases - unsigned int ref_case = (isotropic_only) - ? RefinementCase::isotropic_refinement - : RefinementCase::cut_x; - - for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) - task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case, - fe, matrices[ref_case-1], ref_case, threshold); - - task_group.join_all (); - } - - - - template - void - compute_face_embedding_matrices(const FiniteElement &fe, - FullMatrix (&matrices)[GeometryInfo::max_children_per_face], - const unsigned int face_coarse, - const unsigned int face_fine, - const double threshold) - { - Assert(face_coarse==0, ExcNotImplemented()); - Assert(face_fine==0, ExcNotImplemented()); - - const unsigned int nc = GeometryInfo::max_children_per_face; - const unsigned int n = fe.dofs_per_face; - const unsigned int nd = fe.n_components(); - const unsigned int degree = fe.degree; - - const bool normal = fe.conforms(FiniteElementData::Hdiv); - const bool tangential = fe.conforms(FiniteElementData::Hcurl); - - for (unsigned int i=0; i face_c_dofs(n); - std::vector face_f_dofs(n); - { - unsigned int face_dof=0; - for (unsigned int i=0; i::vertices_per_face; ++i) - { - const unsigned int offset_c = GeometryInfo::face_to_cell_vertices(face_coarse, i) - *fe.dofs_per_vertex; - const unsigned int offset_f = GeometryInfo::face_to_cell_vertices(face_fine, i) - *fe.dofs_per_vertex; - for (unsigned int j=0; j::lines_per_face; ++i) - { - const unsigned int offset_c = fe.first_line_index - + GeometryInfo::face_to_cell_lines(face_coarse, i-1) - *fe.dofs_per_line; - const unsigned int offset_f = fe.first_line_index - + GeometryInfo::face_to_cell_lines(face_fine, i-1) - *fe.dofs_per_line; - for (unsigned int j=0; j::quads_per_face; ++i) - { - const unsigned int offset_c = fe.first_quad_index - + face_coarse - *fe.dofs_per_quad; - const unsigned int offset_f = fe.first_quad_index - + face_fine - *fe.dofs_per_quad; - for (unsigned int j=0; j tria; - GridGenerator::hyper_cube (tria, 0, 1); - tria.refine_global(1); - MappingCartesian mapping; - - // Setup quadrature and FEValues - // for a face. We cannot use - // FEFaceValues and - // FESubfaceValues because of - // some nifty handling of - // refinement cases. Guido stops - // disliking and instead starts - // hating the anisotropic implementation - QGauss q_gauss(degree+1); - const Quadrature q_fine = QProjector::project_to_face(q_gauss, face_fine); - const unsigned int nq = q_fine.size(); - - FEValues fine (mapping, fe, q_fine, - update_quadrature_points | update_JxW_values | update_values); - - // We search for the polynomial on - // the small cell, being equal to - // the coarse polynomial in all - // quadrature points. - - // First build the matrix for this - // least squares problem. This - // contains the values of the fine - // cell polynomials in the fine - // cell grid points. - - // This matrix is the same for all - // children. - fine.reinit(tria.begin_active()); - FullMatrix A(nq*nd, n); - for (unsigned int j=0; j H(A); - - Vector v_coarse(nq*nd); - Vector v_fine(n); - - - - for (unsigned int cell_number = 0; cell_number < GeometryInfo::max_children_per_face; - ++cell_number) - { - const Quadrature q_coarse - = QProjector::project_to_subface(q_gauss, face_coarse, cell_number); - FEValues coarse (mapping, fe, q_coarse, update_values); - - typename Triangulation::active_cell_iterator fine_cell - = tria.begin(0)->child(GeometryInfo::child_cell_on_face( - tria.begin(0)->refinement_case(), face_coarse, cell_number)); - fine.reinit(fine_cell); - coarse.reinit(tria.begin(0)); - - FullMatrix &this_matrix = matrices[cell_number]; - - // Compute this once for each - // coarse grid basis function - for (unsigned int i=0; i - void - compute_projection_matrices(const FiniteElement &fe, - std::vector > > &matrices, - const bool isotropic_only) - { - const unsigned int n = fe.dofs_per_cell; - const unsigned int nd = fe.n_components(); - const unsigned int degree = fe.degree; - - // prepare FEValues, quadrature etc on - // coarse cell - QGauss q_fine(degree+1); - const unsigned int nq = q_fine.size(); - - // create mass matrix on coarse cell. - FullMatrix mass(n, n); - { - // set up a triangulation for coarse cell - Triangulation tr; - GridGenerator::hyper_cube (tr, 0, 1); - - FEValues coarse (fe, q_fine, - update_JxW_values | update_values); - - typename Triangulation::cell_iterator coarse_cell - = tr.begin(0); - coarse.reinit (coarse_cell); - - const std::vector &JxW = coarse.get_JxW_values(); - for (unsigned int i=0; i::isotropic_refinement - : RefinementCase::cut_x; - for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) - { - const unsigned int - nc = GeometryInfo::n_children(RefinementCase(ref_case)); - - for (unsigned int i=0; i tr; - GridGenerator::hyper_cube (tr, 0, 1); - tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); - tr.execute_coarsening_and_refinement(); - - FEValues fine (StaticMappingQ1::mapping, fe, q_fine, - update_quadrature_points | update_JxW_values | - update_values); - - typename Triangulation::cell_iterator coarse_cell - = tr.begin(0); - - Vector v_coarse(n); - Vector v_fine(n); - - for (unsigned int cell_number=0; cell_number &this_matrix = matrices[ref_case-1][cell_number]; - - // Compute right hand side, - // which is a fine level basis - // function tested with the - // coarse level functions. - fine.reinit(coarse_cell->child(cell_number)); - const std::vector > &q_points_fine = fine.get_quadrature_points(); - std::vector > q_points_coarse(q_points_fine.size()); - for (unsigned int q=0; q q_coarse (q_points_coarse, - fine.get_JxW_values()); - FEValues coarse (StaticMappingQ1::mapping, fe, q_coarse, update_values); - coarse.reinit(coarse_cell); - - // Build RHS - - const std::vector &JxW = fine.get_JxW_values(); - - // Outer loop over all fine - // grid shape functions phi_j - for (unsigned int j=0; j - void - add_fe_name(const std::string ¶meter_name, - const FEFactoryBase *factory) - { - // Erase everything after the - // actual class name - std::string name = parameter_name; - unsigned int name_end = - name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); - if (name_end < name.size()) - name.erase(name_end); - // first make sure that no other - // thread intercepts the - // operation of this function; - // for this, acquire the lock - // until we quit this function - Threads::Mutex::ScopedLock lock(fe_name_map_lock); - - Assert(fe_name_map[dim][spacedim].find(name) == fe_name_map[dim][spacedim].end(), - ExcMessage("Cannot change existing element in finite element name list")); - - // Insert the normalized name into - // the map - fe_name_map[dim][spacedim][name] = - std::shared_ptr (factory); - } - - - namespace internal - { - namespace - { - // TODO: this encapsulates the call to the - // dimension-dependent fe_name_map so that we - // have a unique interface. could be done - // smarter? - template - FiniteElement * - get_fe_by_name_ext (std::string &name, - const std::map > - &fe_name_map) - { - // Extract the name of the - // finite element class, which only - // contains characters, numbers and - // underscores. - unsigned int name_end = - name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); - const std::string name_part(name, 0, name_end); - name.erase(0, name_part.size()); - - // now things get a little more - // complicated: FESystem. it's - // more complicated, since we - // have to figure out what the - // base elements are. this can - // only be done recursively - if (name_part == "FESystem") - { - // next we have to get at the - // base elements. start with - // the first. wrap the whole - // block into try-catch to - // make sure we destroy the - // pointers we got from - // recursive calls if one of - // these calls should throw - // an exception - std::vector*> base_fes; - std::vector base_multiplicities; - try - { - // Now, just the [...] - // part should be left. - if (name.size() == 0 || name[0] != '[') - throw (std::string("Invalid first character in ") + name); - do - { - // Erase the - // leading '[' or '-' - name.erase(0,1); - // Now, the name of the - // first base element is - // first... Let's get it - base_fes.push_back (get_fe_by_name_ext (name, - fe_name_map)); - // next check whether - // FESystem placed a - // multiplicity after - // the element name - if (name[0] == '^') - { - // yes. Delete the '^' - // and read this - // multiplicity - name.erase(0,1); - - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - name.erase(0, tmp.second); - // add to length, - // including the '^' - base_multiplicities.push_back (tmp.first); - } - else - // no, so - // multiplicity is - // 1 - base_multiplicities.push_back (1); - - // so that's it for - // this base - // element. base - // elements are - // separated by '-', - // and the list is - // terminated by ']', - // so loop while the - // next character is - // '-' - } - while (name[0] == '-'); - - // so we got to the end - // of the '-' separated - // list. make sure that - // we actually had a ']' - // there - if (name.size() == 0 || name[0] != ']') - throw (std::string("Invalid first character in ") + name); - name.erase(0,1); - // just one more sanity check - Assert ((base_fes.size() == base_multiplicities.size()) - && - (base_fes.size() > 0), - ExcInternalError()); - - // ok, apparently - // everything went ok. so - // generate the composed - // element - FiniteElement *system_element = nullptr; - - // uses new FESystem constructor - // which is independent of - // the number of FEs in the system - system_element = new FESystem(base_fes, base_multiplicities); - - // now we don't need the - // list of base elements - // any more - for (unsigned int i=0; isecond.get(); - const FEFactoryBase *fef=dynamic_cast*>(ptr); - return fef->get(1); - } - else - { - // Make sure no other thread - // is just adding an element - Threads::Mutex::ScopedLock lock (fe_name_map_lock); - AssertThrow (fe_name_map.find(name_part) != fe_name_map.end(), - ExcInvalidFEName(name)); - - // Now, just the (degree) - // or (Quadrature<1>(degree+1)) - // part should be left. - if (name.size() == 0 || name[0] != '(') - throw (std::string("Invalid first character in ") + name); - name.erase(0,1); - if (name[0] != 'Q') - { - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - name.erase(0, tmp.second+1); - const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); - const FEFactoryBase *fef=dynamic_cast*>(ptr); - return fef->get(tmp.first); - } - else - { - unsigned int position = name.find('('); - const std::string quadrature_name(name, 0, position); - name.erase(0,position+1); - if (quadrature_name.compare("QGaussLobatto") == 0) - { - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - // delete "))" - name.erase(0, tmp.second+2); - const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); - const FEFactoryBase *fef=dynamic_cast*>(ptr); - return fef->get(QGaussLobatto<1>(tmp.first)); - } - else if (quadrature_name.compare("QGauss") == 0) - { - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - // delete "))" - name.erase(0, tmp.second+2); - const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); - const FEFactoryBase *fef=dynamic_cast*>(ptr); - return fef->get(QGauss<1>(tmp.first)); - } - else if (quadrature_name.compare("QIterated") == 0) - { - // find sub-quadrature - position = name.find('('); - const std::string subquadrature_name(name, 0, position); - AssertThrow(subquadrature_name.compare("QTrapez") == 0, - ExcNotImplemented("Could not detect quadrature of name " + subquadrature_name)); - // delete "QTrapez()," - name.erase(0,position+3); - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - // delete "))" - name.erase(0, tmp.second+2); - const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); - const FEFactoryBase *fef=dynamic_cast*>(ptr); - return fef->get(QIterated<1>(QTrapez<1>(),tmp.first)); - } - else - { - AssertThrow (false,ExcNotImplemented()); - } - } - } - - - // hm, if we have come thus far, we - // didn't know what to do with the - // string we got. so do as the docs - // say: raise an exception - AssertThrow (false, ExcInvalidFEName(name)); - - // make some compilers happy that - // do not realize that we can't get - // here after throwing - return nullptr; - } - - - - template - FiniteElement *get_fe_by_name (std::string &name) - { - return get_fe_by_name_ext (name, fe_name_map[dim][spacedim]); - } - } - } - - - - template - FiniteElement * - get_fe_from_name (const std::string ¶meter_name) - { - return get_fe_by_name (parameter_name); - } - - - - template - FiniteElement * - get_fe_by_name (const std::string ¶meter_name) - { - std::string name = Utilities::trim(parameter_name); - std::size_t index = 1; - // remove spaces that are not between two word (things that match the - // regular expression [A-Za-z0-9_]) characters. - while (2 < name.size() && index < name.size() - 1) - { - if (name[index] == ' ' && - (!(std::isalnum(name[index - 1]) || name[index - 1] == '_') || - !(std::isalnum(name[index + 1]) || name[index + 1] == '_'))) - { - name.erase(index, 1); - } - else - { - ++index; - } - } - - // Create a version of the name - // string where all template - // parameters are eliminated. - for (unsigned int pos1 = name.find('<'); - pos1 < name.size(); - pos1 = name.find('<')) - { - - const unsigned int pos2 = name.find('>'); - // If there is only a single - // character between those two, - // it should be 'd' or the number - // representing the dimension. - if (pos2-pos1 == 2) - { - const char dimchar = '0' + dim; - (void)dimchar; - if (name.at(pos1+1) != 'd') - Assert (name.at(pos1+1) == dimchar, - ExcInvalidFEDimension(name.at(pos1+1), dim)); - } - else - Assert(pos2-pos1 == 4, ExcInvalidFEName(name)); - - // If pos1==pos2, then we are - // probably at the end of the - // string - if (pos2 != pos1) - name.erase(pos1, pos2-pos1+1); - } - // Replace all occurrences of "^dim" - // by "^d" to be handled by the - // next loop - for (unsigned int pos = name.find("^dim"); - pos < name.size(); - pos = name.find("^dim")) - name.erase(pos+2, 2); - - // Replace all occurrences of "^d" - // by using the actual dimension - for (unsigned int pos = name.find("^d"); - pos < name.size(); - pos = name.find("^d")) - name.at(pos+1) = '0' + dim; - - try - { - FiniteElement *fe = internal::get_fe_by_name (name); - - // Make sure the auxiliary function - // ate up all characters of the name. - AssertThrow (name.size() == 0, - ExcInvalidFEName(parameter_name - + std::string(" extra characters after " - "end of name"))); - return fe; - } - catch (const std::string &errline) - { - AssertThrow(false, ExcInvalidFEName(parameter_name - + std::string(" at ") - + errline)); - return nullptr; - } - } - - - - template - void - compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - FullMatrix &X) - { - Assert (fe.n_components() == 1, ExcNotImplemented()); - - // first build the matrices M and Q - // described in the documentation - FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); - FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); - - for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); - M_inverse.invert (M); - - // finally compute the result - X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); - M_inverse.mmult (X, Q); - - Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); - Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); - } - - - - template - void - compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &quadrature, - FullMatrix &I_q) - { - Assert (fe.n_components() == 1, ExcNotImplemented()); - Assert (I_q.m() == quadrature.size(), - ExcMessage ("Wrong matrix size")); - Assert (I_q.n() == fe.dofs_per_cell, ExcMessage ("Wrong matrix size")); - - for (unsigned int q=0; q - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, - std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes) - { - - // check that the number columns of the projection_matrix - // matches the size of the vector_of_tensors_at_qp - Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), - ExcDimensionMismatch(projection_matrix.n_cols(), - vector_of_tensors_at_qp.size())); - - // check that the number rows of the projection_matrix - // matches the size of the vector_of_tensors_at_nodes - Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), - ExcDimensionMismatch(projection_matrix.n_rows(), - vector_of_tensors_at_nodes.size())); - - // number of support points (nodes) to project to - const unsigned int n_support_points = projection_matrix.n_rows(); - // number of quadrature points to project from - const unsigned int n_quad_points = projection_matrix.n_cols(); - - // component projected to the nodes - Vector component_at_node(n_support_points); - // component at the quadrature point - Vector component_at_qp(n_quad_points); - - for (unsigned int ii = 0; ii < dim; ++ii) - { - - component_at_qp = 0; - - // populate the vector of components at the qps - // from vector_of_tensors_at_qp - // vector_of_tensors_at_qp data is in form: - // columns: 0, 1, ..., dim - // rows: 0,1,...., n_quad_points - // so extract the ii'th column of vector_of_tensors_at_qp - for (unsigned int q = 0; q < n_quad_points; ++q) - { - component_at_qp(q) = vector_of_tensors_at_qp[q][ii]; - } - - // project from the qps -> nodes - // component_at_node = projection_matrix_u * component_at_qp - projection_matrix.vmult(component_at_node, component_at_qp); - - // rewrite the projection of the components - // back into the vector of tensors - for (unsigned int nn =0; nn - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, - std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes) - { - - // check that the number columns of the projection_matrix - // matches the size of the vector_of_tensors_at_qp - Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), - ExcDimensionMismatch(projection_matrix.n_cols(), - vector_of_tensors_at_qp.size())); - - // check that the number rows of the projection_matrix - // matches the size of the vector_of_tensors_at_nodes - Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), - ExcDimensionMismatch(projection_matrix.n_rows(), - vector_of_tensors_at_nodes.size())); - - // number of support points (nodes) - const unsigned int n_support_points = projection_matrix.n_rows(); - // number of quadrature points to project from - const unsigned int n_quad_points = projection_matrix.n_cols(); - - // number of unique entries in a symmetric second-order tensor - const unsigned int n_independent_components = - SymmetricTensor<2, dim >::n_independent_components; - - // component projected to the nodes - Vector component_at_node(n_support_points); - // component at the quadrature point - Vector component_at_qp(n_quad_points); - - // loop over the number of unique dimensions of the tensor - for (unsigned int ii = 0; ii < n_independent_components; ++ii) - { - - component_at_qp = 0; - - // row-column entry of tensor corresponding the unrolled index - TableIndices<2> row_column_index = SymmetricTensor< 2, dim >::unrolled_to_component_indices(ii); - const unsigned int row = row_column_index[0]; - const unsigned int column = row_column_index[1]; - - // populate the vector of components at the qps - // from vector_of_tensors_at_qp - // vector_of_tensors_at_qp is in form: - // columns: 0, 1, ..., n_independent_components - // rows: 0,1,...., n_quad_points - // so extract the ii'th column of vector_of_tensors_at_qp - for (unsigned int q = 0; q < n_quad_points; ++q) - { - component_at_qp(q) = (vector_of_tensors_at_qp[q])[row][column]; - } - - // project from the qps -> nodes - // component_at_node = projection_matrix_u * component_at_qp - projection_matrix.vmult(component_at_node, component_at_qp); - - // rewrite the projection of the components back into the vector of tensors - for (unsigned int nn =0; nn - void - compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face, - FullMatrix &X) - { - Assert (fe.n_components() == 1, ExcNotImplemented()); - Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree)); - - - - // build the matrices M and Q - // described in the documentation - FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); - FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); - - { - // need an FEFaceValues object to evaluate shape function - // values on the specified face. - FEFaceValues fe_face_values (fe, lhs_quadrature, update_values); - fe_face_values.reinit (cell, face); // setup shape_value on this face. - - for (unsigned int i=0; i fe_face_values (fe, rhs_quadrature, update_values); - fe_face_values.reinit (cell, face); // setup shape_value on this face. - - for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); - M_inverse.invert (M); - - // finally compute the result - X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); - M_inverse.mmult (X, Q); - - Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); - Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); - } - - - - template - void - hierarchic_to_lexicographic_numbering (unsigned int degree, std::vector &h2l) - { - // number of support points in each direction - const unsigned int n = degree+1; - - const unsigned int dofs_per_cell = Utilities::fixed_power(n); - - // Assert size maches degree - AssertDimension (h2l.size(), dofs_per_cell); - - // polynomial degree - const unsigned int dofs_per_line = degree - 1; - - // the following lines of code are somewhat odd, due to the way the - // hierarchic numbering is organized. if someone would really want to - // understand these lines, you better draw some pictures where you - // indicate the indices and orders of vertices, lines, etc, along with the - // numbers of the degrees of freedom in hierarchical and lexicographical - // order - switch (dim) - { - case 1: - { - h2l[0] = 0; - h2l[1] = dofs_per_cell-1; - for (unsigned int i=2; i - void - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, - std::vector &h2l) - { - Assert (h2l.size() == fe.dofs_per_cell, - ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); - hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); - } - - - - template - std::vector - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) - { - Assert (fe.n_components() == 1, ExcInvalidFE()); - std::vector h2l(fe.dofs_per_cell); - hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); - return (h2l); - } - - - - template - void - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, - std::vector &l2h) - { - l2h = lexicographic_to_hierarchic_numbering (fe); - } - - - - template - std::vector - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) - { - return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); - } - -} // end of namespace FETools - - - /*-------------- Explicit Instantiations -------------------------------*/ #include "fe_tools.inst" - -/*---------------------------- fe_tools.cc ---------------------------*/ - DEAL_II_NAMESPACE_CLOSE -- 2.39.5