--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/householder.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_q_bubbles.h>
+#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/fe/fe_q_hierarchical.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgp_monomial.h>
+#include <deal.II/fe/fe_dgp_nonparametric.h>
+#include <deal.II/fe/fe_dg_vector.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_abf.h>
+#include <deal.II/fe/fe_bdm.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_rannacher_turek.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_cartesian.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/hp/dof_handler.h>
+
+
+#include <deal.II/base/index_set.h>
+
+#include <cctype>
+#include <iostream>
+#include <memory>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace FETools
+{
+ namespace Compositing
+ {
+
+ template <int dim, int spacedim>
+ FiniteElementData<dim>
+ multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &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; i<fes.size(); i++)
+ if (multiplicities[i]>0)
+ {
+ n_components = fes[i]->n_components();
+ break;
+ }
+
+ for (unsigned int i=0; i<fes.size(); i++)
+ if (multiplicities[i]>0)
+ {
+ 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<dim>::Conformity total_conformity
+ = typename FiniteElementData<dim>::Conformity();
+ {
+ unsigned int index = 0;
+ for (index=0; index<fes.size(); ++index)
+ if (multiplicities[index]>0)
+ {
+ total_conformity = fes[index]->conforming_space;
+ break;
+ }
+
+ for (; index<fes.size(); ++index)
+ if (multiplicities[index]>0)
+ total_conformity =
+ typename FiniteElementData<dim>::Conformity(total_conformity
+ &
+ fes[index]->conforming_space);
+ }
+
+ std::vector<unsigned int> 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<dim> (dpo,
+ (do_tensor_product ? multiplied_n_components : n_components),
+ degree,
+ total_conformity,
+ block_indices);
+ }
+
+
+
+ template <int dim, int spacedim>
+ FiniteElementData<dim>
+ multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2,
+ const unsigned int N2,
+ const FiniteElement<dim,spacedim> *fe3,
+ const unsigned int N3,
+ const FiniteElement<dim,spacedim> *fe4,
+ const unsigned int N4,
+ const FiniteElement<dim,spacedim> *fe5,
+ const unsigned int N5)
+ {
+ std::vector<const FiniteElement<dim,spacedim>*> fes;
+ fes.push_back(fe1);
+ fes.push_back(fe2);
+ fes.push_back(fe3);
+ fes.push_back(fe4);
+ fes.push_back(fe5);
+
+ std::vector<unsigned int> 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 <int dim, int spacedim>
+ std::vector<bool>
+ compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) // 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<bool> 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<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ 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 (index_in_base < fes[base]->dofs_per_cell,
+ ExcInternalError());
+ retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+ }
+ }
+
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ 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 (index_in_base < fes[base]->dofs_per_cell,
+ ExcInternalError());
+ retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+ }
+ }
+
+ // 3. Quads
+ if (GeometryInfo<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ 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 (index_in_base < fes[base]->dofs_per_cell,
+ ExcInternalError());
+ retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+ }
+ }
+
+ // 4. Hexes
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ 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 (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 <int dim, int spacedim>
+ std::vector<bool>
+ compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2,
+ const unsigned int N2,
+ const FiniteElement<dim,spacedim> *fe3,
+ const unsigned int N3,
+ const FiniteElement<dim,spacedim> *fe4,
+ const unsigned int N4,
+ const FiniteElement<dim,spacedim> *fe5,
+ const unsigned int N5)
+ {
+ std::vector<const FiniteElement<dim,spacedim>*> fe_list;
+ std::vector<unsigned int> 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 <int dim, int spacedim>
+ std::vector<ComponentMask>
+ compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
+ const std::vector<unsigned int> &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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //needed because fe might be NULL
+ n_components += fes[i]->n_components() * multiplicities[i];
+ }
+ else
+ {
+ for (unsigned int i=0; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //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; i<fes.size(); ++i)
+ if (multiplicities[i]>0) //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<std::vector<bool> >
+ retval (n_shape_functions, std::vector<bool> (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<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<vector<bool> > into a vector<ComponentMask>.
+ // 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<bool> anyway.
+ std::vector<ComponentMask> xretval (retval.size());
+ for (unsigned int i=0; i<retval.size(); ++i)
+ xretval[i] = ComponentMask(retval[i]);
+ return xretval;
+ }
+
+
+
+ /**
+ * Compute the non-zero vector components of a composed finite element.
+ */
+ template <int dim, int spacedim>
+ std::vector<ComponentMask>
+ compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
+ const unsigned int N1,
+ const FiniteElement<dim,spacedim> *fe2,
+ const unsigned int N2,
+ const FiniteElement<dim,spacedim> *fe3,
+ const unsigned int N3,
+ const FiniteElement<dim,spacedim> *fe4,
+ const unsigned int N4,
+ const FiniteElement<dim,spacedim> *fe5,
+ const unsigned int N5,
+ const bool do_tensor_product)
+ {
+ std::vector<const FiniteElement<dim,spacedim>*> fe_list;
+ std::vector<unsigned int> 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 <int dim, int spacedim>
+ 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<dim,spacedim> &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<fe.base_element(base).n_components(); ++k)
+ component_to_base_table[total_index++]
+ = std::make_pair(std::make_pair(base,k), m);
+ }
+ Assert (total_index == component_to_base_table.size(),
+ ExcInternalError());
+ }
+ else
+ {
+ // The base element establishing a component does not make sense in this case.
+ // Set up to something meaningless:
+ for (unsigned int i = 0; i < component_to_base_table.size(); i++)
+ component_to_base_table[i] = std::make_pair(std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int), numbers::invalid_unsigned_int);
+
+ }
+
+
+ // Initialize index tables. Multi-component base elements have to be
+ // thought of. For non-primitive shape functions, have a special invalid
+ // index.
+ const std::pair<unsigned int, unsigned int>
+ 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<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_vertex;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base, m), index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = fe.base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
+ }
+ else
+ system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_line;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_line*line_number +
+ local_index +
+ fe.base_element(base).first_line_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = fe.base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
+ }
+ else
+ system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+
+ // 3. Quads
+ if (GeometryInfo<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_quad;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_quad*quad_number +
+ local_index +
+ fe.base_element(base).first_quad_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = fe.base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
+ }
+ else
+ system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+
+ // 4. Hexes
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_hex;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_hex*hex_number +
+ local_index +
+ fe.base_element(base).first_hex_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = fe.base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
+ }
+ else
+ system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ 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<dim,spacedim> &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<unsigned int, unsigned int>
+ 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<GeometryInfo<dim>::vertices_per_face;
+ ++vertex_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_vertex;
+ ++local_index, ++total_index)
+ {
+ // get (cell) index of this shape function inside the base
+ // element to see whether the shape function is primitive
+ // (assume that all shape functions on vertices share the same
+ // primitivity property; assume likewise for all shape functions
+ // located on lines, quads, etc. this way, we can ask for
+ // primitivity of only _one_ shape function, which is taken as
+ // representative for all others located on the same type of
+ // object):
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
+
+ const unsigned int face_index_in_base
+ = (fe.base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
+ }
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_face > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::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<fe.element_multiplicity(base);
+ ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_line;
+ ++local_index, ++total_index)
+ {
+ // do everything alike for this type of object
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_line*line_number +
+ local_index +
+ fe.base_element(base).first_line_index);
+
+ const unsigned int face_index_in_base
+ = (fe.base_element(base).first_face_line_index +
+ fe.base_element(base).dofs_per_line * line_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
+ }
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+
+ // 3. Quads
+ if (GeometryInfo<dim>::quads_per_face > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_face;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fe.n_base_elements(); ++base)
+ for (unsigned int m=0; m<fe.element_multiplicity(base);
+ ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
+ for (unsigned int local_index = 0;
+ local_index < fe.base_element(base).dofs_per_quad;
+ ++local_index, ++total_index)
+ {
+ // do everything alike for this type of object
+ const unsigned int index_in_base
+ = (fe.base_element(base).dofs_per_quad*quad_number +
+ local_index +
+ fe.base_element(base).first_quad_index);
+
+ const unsigned int face_index_in_base
+ = (fe.base_element(base).first_face_quad_index +
+ fe.base_element(base).dofs_per_quad * quad_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
+ }
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+ Assert (total_index == fe.dofs_per_face, ExcInternalError());
+ Assert (total_index == face_system_to_component_table.size(),
+ ExcInternalError());
+ Assert (total_index == face_system_to_base_table.size(),
+ ExcInternalError());
+ }
+ }
+
+
+
+ // Not implemented in the general case.
+ template <class FE>
+ FiniteElement<FE::dimension, FE::space_dimension> *
+ FEFactory<FE>::get (const Quadrature<1> &) const
+ {
+ Assert(false, ExcNotImplemented());
+ return nullptr;
+ }
+
+ // Specializations for FE_Q.
+ template <>
+ FiniteElement<1, 1> *
+ FEFactory<FE_Q<1, 1> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q<1>(quad);
+ }
+
+ template <>
+ FiniteElement<2, 2> *
+ FEFactory<FE_Q<2, 2> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q<2>(quad);
+ }
+
+ template <>
+ FiniteElement<3, 3> *
+ FEFactory<FE_Q<3, 3> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q<3>(quad);
+ }
+
+ // Specializations for FE_Q_DG0.
+ template <>
+ FiniteElement<1, 1> *
+ FEFactory<FE_Q_DG0<1, 1> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q_DG0<1>(quad);
+ }
+
+ template <>
+ FiniteElement<2, 2> *
+ FEFactory<FE_Q_DG0<2, 2> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q_DG0<2>(quad);
+ }
+
+ template <>
+ FiniteElement<3, 3> *
+ FEFactory<FE_Q_DG0<3, 3> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q_DG0<3>(quad);
+ }
+
+ // Specializations for FE_Q_Bubbles.
+ template <>
+ FiniteElement<1, 1> *
+ FEFactory<FE_Q_Bubbles<1, 1> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q_Bubbles<1>(quad);
+ }
+
+ template <>
+ FiniteElement<2, 2> *
+ FEFactory<FE_Q_Bubbles<2, 2> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q_Bubbles<2>(quad);
+ }
+
+ template <>
+ FiniteElement<3, 3> *
+ FEFactory<FE_Q_Bubbles<3, 3> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_Q_Bubbles<3>(quad);
+ }
+
+ // Specializations for FE_DGQArbitraryNodes.
+ template <>
+ FiniteElement<1, 1> *
+ FEFactory<FE_DGQ<1> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_DGQArbitraryNodes<1>(quad);
+ }
+
+ template <>
+ FiniteElement<1, 2> *
+ FEFactory<FE_DGQ<1, 2> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_DGQArbitraryNodes<1, 2>(quad);
+ }
+
+ template <>
+ FiniteElement<1, 3> *
+ FEFactory<FE_DGQ<1, 3> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_DGQArbitraryNodes<1, 3>(quad);
+ }
+
+ template <>
+ FiniteElement<2, 2> *
+ FEFactory<FE_DGQ<2> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_DGQArbitraryNodes<2>(quad);
+ }
+
+ template <>
+ FiniteElement<2, 3> *
+ FEFactory<FE_DGQ<2, 3> >::get (const Quadrature<1> &quad) const
+ {
+ return new FE_DGQArbitraryNodes<2, 3>(quad);
+ }
+
+ template <>
+ FiniteElement<3, 3> *
+ FEFactory<FE_DGQ<3> >::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 <int dim>
+ void
+ fill_no_codim_fe_names (std::map<std::string,std::shared_ptr<const Subscriptor> > &result)
+ {
+ typedef std::shared_ptr<const Subscriptor> FEFactoryPointer;
+
+ result["FE_Q_Hierarchical"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q_Hierarchical<dim> >);
+ result["FE_ABF"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_ABF<dim> >);
+ result["FE_Bernstein"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Bernstein<dim> >);
+ result["FE_BDM"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_BDM<dim> >);
+ result["FE_DGBDM"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGBDM<dim> >);
+ result["FE_DGNedelec"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGNedelec<dim> >);
+ result["FE_DGRaviartThomas"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGRaviartThomas<dim> >);
+ result["FE_RaviartThomas"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomas<dim> >);
+ result["FE_RaviartThomasNodal"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomasNodal<dim> >);
+ result["FE_Nedelec"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Nedelec<dim> >);
+ result["FE_DGPNonparametric"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGPNonparametric<dim> >);
+ result["FE_DGP"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGP<dim> >);
+ result["FE_DGPMonomial"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGPMonomial<dim> >);
+ result["FE_DGQ"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim> >);
+ result["FE_DGQArbitraryNodes"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim> >);
+ result["FE_DGQLegendre"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQLegendre<dim> >);
+ result["FE_DGQHermite"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQHermite<dim> >);
+ result["FE_FaceQ"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_FaceQ<dim> >);
+ result["FE_FaceP"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_FaceP<dim> >);
+ result["FE_Q"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q<dim> >);
+ result["FE_Q_DG0"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q_DG0<dim> >);
+ result["FE_Q_Bubbles"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q_Bubbles<dim> >);
+ result["FE_Q_iso_Q1"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q_iso_Q1<dim> >);
+ result["FE_Nothing"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Nothing<dim> >);
+ result["FE_RannacherTurek"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_RannacherTurek<dim> >);
+ }
+
+
+
+ // This function fills a map from names to finite elements for any
+ // dimension and codimension for those elements which support
+ // nonzero codimension.
+ template <int dim, int spacedim>
+ void
+ fill_codim_fe_names (std::map<std::string,std::shared_ptr<const Subscriptor> > &result)
+ {
+ typedef std::shared_ptr<const Subscriptor> FEFactoryPointer;
+
+ result["FE_Bernstein"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Bernstein<dim,spacedim> >);
+ result["FE_DGP"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGP<dim,spacedim> >);
+ result["FE_DGQ"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim,spacedim> >);
+ result["FE_Nothing"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Nothing<dim,spacedim> >);
+ result["FE_DGQArbitraryNodes"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim,spacedim> >);
+ result["FE_DGQLegendre"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQLegendre<dim,spacedim> >);
+ result["FE_DGQHermite"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_DGQHermite<dim,spacedim> >);
+ result["FE_Q_Bubbles"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q_Bubbles<dim,spacedim> >);
+ result["FE_Q_DG0"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q_DG0<dim,spacedim> >);
+ result["FE_Q_iso_Q1"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q_iso_Q1<dim,spacedim> >);
+ result["FE_Q"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Q<dim,spacedim> >);
+ result["FE_Bernstein"]
+ = FEFactoryPointer(new FETools::FEFactory<FE_Bernstein<dim,spacedim> >);
+ }
+
+ // 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<std::vector<
+ std::map<std::string,
+ std::shared_ptr<const Subscriptor> > > >
+ fill_default_map()
+ {
+ std::vector<std::vector<
+ std::map<std::string,
+ std::shared_ptr<const Subscriptor> > > >
+ 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<dim,spacedim> 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<std::vector<
+ std::map<std::string,
+ std::shared_ptr<const Subscriptor> > > >
+ fe_name_map = fill_default_map();
+}
+
+
+
+
+
+
+namespace
+{
+
+ // forwarder function for
+ // FE::get_interpolation_matrix. we
+ // will want to call that function
+ // for arbitrary FullMatrix<T>
+ // 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 <int dim, int spacedim>
+ inline
+ void gim_forwarder (const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<double> &interpolation_matrix)
+ {
+ fe2.get_interpolation_matrix (fe1, interpolation_matrix);
+ }
+
+
+
+ template <int dim, typename number, int spacedim>
+ inline
+ void gim_forwarder (const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &interpolation_matrix)
+ {
+ FullMatrix<double> 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 "<dim>" or the
+ // specialized string with "dim"
+ // replaced with the numeric value
+ // of the template argument
+ template <int dim, int spacedim>
+ 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 <int dim, int spacedim>
+ FEFactoryBase<dim,spacedim>::~FEFactoryBase()
+ {}
+
+
+
+ template <int dim, int spacedim>
+ void compute_component_wise(
+ const FiniteElement<dim,spacedim> &element,
+ std::vector<unsigned int> &renumbering,
+ std::vector<std::vector<unsigned int> > &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<comp_start.size(); ++i)
+ {
+ comp_start[i].resize(element.element_multiplicity(i));
+ const unsigned int increment
+ = element.base_element(i).dofs_per_cell;
+
+ for (unsigned int j=0; j<comp_start[i].size(); ++j)
+ {
+ comp_start[i][j] = k;
+ k += increment;
+ }
+ }
+
+ // For each index i of the
+ // unstructured cellwise
+ // numbering, renumbering
+ // contains the index of the
+ // cell-block numbering
+ for (unsigned int i=0; i<element.dofs_per_cell; ++i)
+ {
+ std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
+ indices = element.system_to_base_index(i);
+ renumbering[i] = comp_start[indices.first.first][indices.first.second]
+ +indices.second;
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ void compute_block_renumbering (
+ const FiniteElement<dim,spacedim> &element,
+ std::vector<types::global_dof_index> &renumbering,
+ std::vector<types::global_dof_index> &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<element.n_base_elements(); ++b)
+ for (unsigned int m=0; m<element.element_multiplicity(b); ++m)
+ {
+ block_data[count++] = (return_start_indices)
+ ? k
+ : (element.base_element(b).n_dofs_per_cell());
+ k += element.base_element(b).n_dofs_per_cell();
+ }
+ Assert (count == element.n_blocks(), ExcInternalError());
+
+ std::vector<types::global_dof_index> start_indices(block_data.size());
+ k = 0;
+ for (unsigned int i=0; i<block_data.size(); ++i)
+ if (return_start_indices)
+ start_indices[i] = block_data[i];
+ else
+ {
+ start_indices[i] = k;
+ k += block_data[i];
+ }
+
+ for (unsigned int i=0; i<element.dofs_per_cell; ++i)
+ {
+ std::pair<unsigned int, types::global_dof_index>
+ indices = element.system_to_block_index(i);
+ renumbering[i] = start_indices[indices.first]
+ +indices.second;
+ }
+ }
+
+
+
+ template <int dim, typename number, int spacedim>
+ void get_interpolation_matrix (const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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<dim,spacedim>::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<Point<dim> > &
+ fe2_support_points = fe2.get_unit_support_points ();
+
+ typedef FiniteElement<dim,spacedim> FEL;
+ Assert(fe2_support_points.size()==fe2.dofs_per_cell,
+ typename FEL::ExcFEHasNoSupportPoints());
+
+ for (unsigned int i=0; i<fe2.dofs_per_cell; ++i)
+ {
+ const unsigned int i1 = fe2.system_to_component_index(i).first;
+ for (unsigned int j=0; j<fe1.dofs_per_cell; ++j)
+ {
+ const unsigned int j1 = fe1.system_to_component_index(j).first;
+ if (i1==j1)
+ interpolation_matrix(i,j) = fe1.shape_value (j,fe2_support_points[i]);
+ else
+ interpolation_matrix(i,j)=0.;
+ }
+ }
+ }
+
+
+
+ template <int dim, typename number, int spacedim>
+ void get_back_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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<number> first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell);
+ FullMatrix<number> 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 <int dim, typename number, int spacedim>
+ void get_interpolation_difference_matrix (const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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<number> interpolation_matrix(fe1.dofs_per_cell);
+ get_back_interpolation_matrix(fe1, fe2, interpolation_matrix);
+
+ for (unsigned int i=0; i<fe1.dofs_per_cell; ++i)
+ difference_matrix(i,i) = 1.;
+
+ // compute difference
+ difference_matrix.add (-1, interpolation_matrix);
+ }
+
+
+
+ template <int dim, typename number, int spacedim>
+ void get_projection_matrix (const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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<dim,spacedim> 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<dim> quadrature(degree+1);
+ // Set up FEValues.
+ const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values;
+ FEValues<dim> val1 (fe1, quadrature, update_values);
+ val1.reinit (tr.begin_active());
+ FEValues<dim> val2 (fe2, quadrature, flags);
+ val2.reinit (tr.begin_active());
+
+ // Integrate and invert mass matrix
+ // This happens in the target space
+ FullMatrix<double> mass (n2, n2);
+
+ for (unsigned int k=0; k<quadrature.size(); ++k)
+ {
+ const double w = val2.JxW(k);
+ for (unsigned int i=0; i<n2; ++i)
+ {
+ const double v = val2.shape_value(i,k);
+ for (unsigned int j=0; j<n2; ++j)
+ mass(i,j) += w*v * val2.shape_value(j,k);
+ }
+ }
+ // Gauss-Jordan should be
+ // sufficient since we expect the
+ // mass matrix to be
+ // well-conditioned
+ mass.gauss_jordan();
+
+ // Now, test every function of fe1
+ // with test functions of fe2 and
+ // compute the projection of each
+ // unit vector.
+ Vector<double> b(n2);
+ Vector<double> x(n2);
+
+ for (unsigned int j=0; j<n1; ++j)
+ {
+ b = 0.;
+ for (unsigned int i=0; i<n2; ++i)
+ for (unsigned int k=0; k<quadrature.size(); ++k)
+ {
+ const double w = val2.JxW(k);
+ const double u = val1.shape_value(j,k);
+ const double v = val2.shape_value(i,k);
+ b(i) += u*v*w;
+ }
+
+ // Multiply by the inverse
+ mass.vmult(x,b);
+ for (unsigned int i=0; i<n2; ++i)
+ matrix(i,j) = x(i);
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ FullMatrix<double>
+ compute_node_matrix(const FiniteElement<dim,spacedim> &fe)
+ {
+ const unsigned int n_dofs = fe.dofs_per_cell;
+
+ FullMatrix<double> N (n_dofs, n_dofs);
+
+ Assert (fe.has_generalized_support_points(), ExcNotInitialized());
+ Assert (fe.n_components() == dim, ExcNotImplemented());
+
+ const std::vector<Point<dim> > &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<Vector<double> >
+ support_point_values (points.size(), Vector<double>(dim));
+
+ // In this vector, we store the
+ // result of the interpolation
+ std::vector<double> 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<n_dofs; ++i)
+ {
+ // get the values of the current set of shape functions
+ // at the generalized support points
+ for (unsigned int k=0; k<points.size(); ++k)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ support_point_values[k][d] = fe.shape_value_component(i, points[k], d);
+ Assert (numbers::is_finite(support_point_values[k][d]), ExcInternalError());
+ }
+
+ fe.convert_generalized_support_point_values_to_dof_values(support_point_values,
+ nodal_values);
+
+ // Enter the interpolated dofs into the matrix
+ for (unsigned int j=0; j<n_dofs; ++j)
+ {
+ N(j,i) = nodal_values[j];
+ Assert (numbers::is_finite(nodal_values[j]), ExcInternalError());
+ }
+ }
+
+ return N;
+ }
+
+
+
+ /*
+ template <>
+ void
+ compute_embedding_matrices(const FiniteElement<1,2> &,
+ std::vector<std::vector<FullMatrix<double> > > &,
+ const bool)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+
+ template <>
+ void
+ compute_embedding_matrices(const FiniteElement<1,3> &,
+ std::vector<std::vector<FullMatrix<double> > > &,
+ const bool)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+
+
+ template <>
+ void
+ compute_embedding_matrices(const FiniteElement<2,3>&,
+ std::vector<std::vector<FullMatrix<double> > >&,
+ const bool)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+ */
+
+ namespace
+ {
+ template <int dim, typename number, int spacedim>
+ void
+ compute_embedding_for_shape_function (
+ const unsigned int i,
+ const FiniteElement<dim, spacedim> &fe,
+ const FEValues<dim, spacedim> &coarse,
+ const Householder<double> &H,
+ FullMatrix<number> &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<number> v_coarse(nq*nd);
+ Vector<number> 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 <int dim, typename number, int spacedim>
+ void
+ compute_embedding_matrices_for_refinement_case (
+ const FiniteElement<dim, spacedim> &fe,
+ std::vector<FullMatrix<number> > &matrices,
+ const unsigned int ref_case,
+ const double threshold)
+ {
+ const unsigned int n = fe.dofs_per_cell;
+ const unsigned int nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(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<dim,spacedim> tria;
+ GridGenerator::hyper_cube (tria, 0, 1);
+ tria.begin_active()->set_refine_flag (RefinementCase<dim>(ref_case));
+ tria.execute_coarsening_and_refinement ();
+
+ const unsigned int degree = fe.degree;
+ QGauss<dim> q_fine (degree+1);
+ const unsigned int nq = q_fine.size();
+
+ FEValues<dim,spacedim> 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<number> 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<double> H (A);
+ unsigned int cell_number = 0;
+
+ Threads::TaskGroup<void> task_group;
+
+ for (typename Triangulation<dim,spacedim>::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<Point<spacedim> > &q_points_fine = fine.get_quadrature_points();
+ std::vector<Point<dim> > q_points_coarse(q_points_fine.size());
+ for (unsigned int i=0; i<q_points_fine.size(); ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ q_points_coarse[i](j) = q_points_fine[i](j);
+ const Quadrature<dim> q_coarse (q_points_coarse,
+ fine.get_JxW_values ());
+ FEValues<dim,spacedim> coarse (fe, q_coarse, update_values);
+
+ coarse.reinit (tria.begin (0));
+
+ FullMatrix<double> &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<dim, number, spacedim>,
+ 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<dim, number, spacedim>
+ (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<dim>::n_children (RefinementCase<dim> (ref_case)),
+ ExcInternalError ());
+ }
+ }
+
+
+
+ template <int dim, typename number, int spacedim>
+ void
+ compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
+ std::vector<std::vector<FullMatrix<number> > > &matrices,
+ const bool isotropic_only,
+ const double threshold)
+ {
+ Threads::TaskGroup<void> task_group;
+
+ // loop over all possible refinement cases
+ unsigned int ref_case = (isotropic_only)
+ ? RefinementCase<dim>::isotropic_refinement
+ : RefinementCase<dim>::cut_x;
+
+ for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
+ task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case<dim, number, spacedim>,
+ fe, matrices[ref_case-1], ref_case, threshold);
+
+ task_group.join_all ();
+ }
+
+
+
+ template <int dim, typename number, int spacedim>
+ void
+ compute_face_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
+ FullMatrix<number> (&matrices)[GeometryInfo<dim>::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<dim>::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<dim>::Hdiv);
+ const bool tangential = fe.conforms(FiniteElementData<dim>::Hcurl);
+
+ 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));
+ }
+
+ // In order to make the loops below
+ // simpler, we introduce vectors
+ // containing for indices 0-n the
+ // number of the corresponding
+ // shape value on the cell.
+ std::vector<unsigned int> face_c_dofs(n);
+ std::vector<unsigned int> face_f_dofs(n);
+ {
+ unsigned int face_dof=0;
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
+ {
+ const unsigned int offset_c = GeometryInfo<dim>::face_to_cell_vertices(face_coarse, i)
+ *fe.dofs_per_vertex;
+ const unsigned int offset_f = GeometryInfo<dim>::face_to_cell_vertices(face_fine, i)
+ *fe.dofs_per_vertex;
+ for (unsigned int j=0; j<fe.dofs_per_vertex; ++j)
+ {
+ face_c_dofs[face_dof] = offset_c + j;
+ face_f_dofs[face_dof] = offset_f + j;
+ ++face_dof;
+ }
+ }
+ for (unsigned int i=1; i<=GeometryInfo<dim>::lines_per_face; ++i)
+ {
+ const unsigned int offset_c = fe.first_line_index
+ + GeometryInfo<dim>::face_to_cell_lines(face_coarse, i-1)
+ *fe.dofs_per_line;
+ const unsigned int offset_f = fe.first_line_index
+ + GeometryInfo<dim>::face_to_cell_lines(face_fine, i-1)
+ *fe.dofs_per_line;
+ for (unsigned int j=0; j<fe.dofs_per_line; ++j)
+ {
+ face_c_dofs[face_dof] = offset_c + j;
+ face_f_dofs[face_dof] = offset_f + j;
+ ++face_dof;
+ }
+ }
+ for (unsigned int i=1; i<=GeometryInfo<dim>::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<fe.dofs_per_quad; ++j)
+ {
+ face_c_dofs[face_dof] = offset_c + j;
+ face_f_dofs[face_dof] = offset_f + j;
+ ++face_dof;
+ }
+ }
+ Assert (face_dof == fe.dofs_per_face, ExcInternalError());
+ }
+
+ // Set up meshes, one with a single
+ // reference cell and refine it once
+ Triangulation<dim,spacedim> tria;
+ GridGenerator::hyper_cube (tria, 0, 1);
+ tria.refine_global(1);
+ MappingCartesian<dim> 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<dim-1> q_gauss(degree+1);
+ const Quadrature<dim> q_fine = QProjector<dim>::project_to_face(q_gauss, face_fine);
+ const unsigned int nq = q_fine.size();
+
+ FEValues<dim> 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<number> A(nq*nd, n);
+ for (unsigned int j=0; j<n; ++j)
+ for (unsigned int k=0; k<nq; ++k)
+ if (nd != dim)
+ for (unsigned int d=0; d<nd; ++d)
+ A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
+ else
+ {
+ if (normal)
+ A(k*nd,j) = fine.shape_value_component(face_f_dofs[j],k,0);
+ if (tangential)
+ for (unsigned int d=1; d<dim; ++d)
+ A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
+ }
+
+ Householder<double> H(A);
+
+ Vector<number> v_coarse(nq*nd);
+ Vector<number> v_fine(n);
+
+
+
+ for (unsigned int cell_number = 0; cell_number < GeometryInfo<dim>::max_children_per_face;
+ ++cell_number)
+ {
+ const Quadrature<dim> q_coarse
+ = QProjector<dim>::project_to_subface(q_gauss, face_coarse, cell_number);
+ FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
+
+ typename Triangulation<dim,spacedim>::active_cell_iterator fine_cell
+ = tria.begin(0)->child(GeometryInfo<dim>::child_cell_on_face(
+ tria.begin(0)->refinement_case(), face_coarse, cell_number));
+ fine.reinit(fine_cell);
+ coarse.reinit(tria.begin(0));
+
+ FullMatrix<double> &this_matrix = matrices[cell_number];
+
+ // Compute this once for each
+ // coarse grid basis function
+ for (unsigned int i=0; i<n; ++i)
+ {
+ // The right hand side of
+ // the least squares
+ // problem consists of the
+ // function values of the
+ // coarse grid function in
+ // each quadrature point.
+ for (unsigned int k=0; k<nq; ++k)
+ if (nd != dim)
+ for (unsigned int d=0; d<nd; ++d)
+ v_coarse(k*nd+d) = coarse.shape_value_component (face_c_dofs[i],k,d);
+ else
+ {
+ if (normal)
+ v_coarse(k*nd) = coarse.shape_value_component(face_c_dofs[i],k,0);
+ if (tangential)
+ for (unsigned int d=1; d<dim; ++d)
+ v_coarse(k*nd+d) = coarse.shape_value_component(face_c_dofs[i],k,d);
+ }
+ // solve the least squares
+ // problem.
+ const double result = H.least_squares(v_fine, v_coarse);
+ Assert (result <= threshold, ExcLeastSquaresError(result));
+ // Avoid compiler 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);
+ }
+ // 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.;
+ }
+ }
+
+
+
+ template <int dim, typename number, int spacedim>
+ void
+ compute_projection_matrices(const FiniteElement<dim,spacedim> &fe,
+ std::vector<std::vector<FullMatrix<number> > > &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<dim> q_fine(degree+1);
+ const unsigned int nq = q_fine.size();
+
+ // create mass matrix on coarse cell.
+ FullMatrix<number> mass(n, n);
+ {
+ // set up a triangulation for coarse cell
+ Triangulation<dim,spacedim> tr;
+ GridGenerator::hyper_cube (tr, 0, 1);
+
+ FEValues<dim,spacedim> coarse (fe, q_fine,
+ update_JxW_values | update_values);
+
+ typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
+ = tr.begin(0);
+ coarse.reinit (coarse_cell);
+
+ const std::vector<double> &JxW = coarse.get_JxW_values();
+ for (unsigned int i=0; i<n; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ if (fe.is_primitive())
+ {
+ const double *coarse_i = &coarse.shape_value(i,0);
+ const double *coarse_j = &coarse.shape_value(j,0);
+ double mass_ij = 0;
+ for (unsigned int k=0; k<nq; ++k)
+ mass_ij += JxW[k] * coarse_i[k] * coarse_j[k];
+ mass(i,j) = mass_ij;
+ }
+ else
+ {
+ double mass_ij = 0;
+ for (unsigned int d=0; d<nd; ++d)
+ for (unsigned int k=0; k<nq; ++k)
+ mass_ij += JxW[k] * coarse.shape_value_component(i,k,d)
+ * coarse.shape_value_component(j,k,d);
+ mass(i,j) = mass_ij;
+ }
+
+ // invert mass matrix
+ mass.gauss_jordan();
+ }
+
+ // loop over all possible
+ // refinement cases
+ unsigned int ref_case = (isotropic_only)
+ ? RefinementCase<dim>::isotropic_refinement
+ : RefinementCase<dim>::cut_x;
+ for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
+ {
+ const unsigned int
+ nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
+
+ for (unsigned int i=0; i<nc; ++i)
+ {
+ Assert(matrices[ref_case-1][i].n() == n,
+ ExcDimensionMismatch(matrices[ref_case-1][i].n(),n));
+ Assert(matrices[ref_case-1][i].m() == n,
+ ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
+ }
+
+ // create a respective refinement on the
+ // triangulation
+ Triangulation<dim,spacedim> tr;
+ GridGenerator::hyper_cube (tr, 0, 1);
+ tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
+ tr.execute_coarsening_and_refinement();
+
+ FEValues<dim,spacedim> fine (StaticMappingQ1<dim,spacedim>::mapping, fe, q_fine,
+ update_quadrature_points | update_JxW_values |
+ update_values);
+
+ typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
+ = tr.begin(0);
+
+ Vector<number> v_coarse(n);
+ Vector<number> v_fine(n);
+
+ for (unsigned int cell_number=0; cell_number<nc; ++cell_number)
+ {
+ FullMatrix<double> &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<Point<spacedim> > &q_points_fine = fine.get_quadrature_points();
+ std::vector<Point<dim> > q_points_coarse(q_points_fine.size());
+ for (unsigned int q=0; q<q_points_fine.size(); ++q)
+ for (unsigned int j=0; j<dim; ++j)
+ q_points_coarse[q](j) = q_points_fine[q](j);
+ Quadrature<dim> q_coarse (q_points_coarse,
+ fine.get_JxW_values());
+ FEValues<dim,spacedim> coarse (StaticMappingQ1<dim,spacedim>::mapping, fe, q_coarse, update_values);
+ coarse.reinit(coarse_cell);
+
+ // Build RHS
+
+ const std::vector<double> &JxW = fine.get_JxW_values();
+
+ // Outer loop over all fine
+ // grid shape functions phi_j
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ {
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ if (fe.is_primitive())
+ {
+ const double *coarse_i = &coarse.shape_value(i,0);
+ const double *fine_j = &fine.shape_value(j,0);
+
+ double update = 0;
+ for (unsigned int k=0; k<nq; ++k)
+ update += JxW[k] * coarse_i[k] * fine_j[k];
+ v_fine(i) = update;
+ }
+ else
+ {
+ double update = 0;
+ for (unsigned int d=0; d<nd; ++d)
+ for (unsigned int k=0; k<nq; ++k)
+ update += JxW[k] * coarse.shape_value_component(i,k,d)
+ * fine.shape_value_component(j,k,d);
+ v_fine(i) = update;
+ }
+ }
+
+ // RHS ready. Solve system
+ // and enter row into
+ // matrix
+ mass.vmult (v_coarse, v_fine);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ this_matrix(i,j) = v_coarse(i);
+ }
+
+ // 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.;
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ add_fe_name(const std::string ¶meter_name,
+ const FEFactoryBase<dim,spacedim> *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<const Subscriptor> (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 <int dim, int spacedim>
+ FiniteElement<dim,spacedim> *
+ get_fe_by_name_ext (std::string &name,
+ const std::map<std::string,
+ std::shared_ptr<const Subscriptor> >
+ &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<const FiniteElement<dim,spacedim>*> base_fes;
+ std::vector<unsigned int> 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<dim,spacedim> (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<int,unsigned int> 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<dim,spacedim> *system_element = nullptr;
+
+ // uses new FESystem constructor
+ // which is independent of
+ // the number of FEs in the system
+ system_element = new FESystem<dim,spacedim>(base_fes, base_multiplicities);
+
+ // now we don't need the
+ // list of base elements
+ // any more
+ for (unsigned int i=0; i<base_fes.size(); ++i)
+ delete base_fes[i];
+
+ // finally return our
+ // findings
+ // Add the closing ']' to
+ // the length
+ return system_element;
+
+ }
+ catch (...)
+ {
+ // ups, some exception
+ // was thrown. prevent a
+ // memory leak, and then
+ // pass on the exception
+ // to the caller
+ for (unsigned int i=0; i<base_fes.size(); ++i)
+ delete base_fes[i];
+ throw;
+ }
+
+ // this is a place where we
+ // should really never get,
+ // since above we have either
+ // returned from the
+ // try-clause, or have
+ // re-thrown in the catch
+ // clause. check that we
+ // never get here
+ Assert (false, ExcInternalError());
+ }
+ else if (name_part == "FE_Nothing")
+ {
+ // remove the () from FE_Nothing()
+ name.erase(0,2);
+
+ // this is a bit of a hack, as
+ // FE_Nothing does not take a
+ // degree, but it does take an
+ // argument, which defaults to 1,
+ // so this properly returns
+ // FE_Nothing()
+ const Subscriptor *ptr = fe_name_map.find(name_part)->second.get();
+ const FEFactoryBase<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(ptr);
+ return fef->get(QGaussLobatto<1>(tmp.first));
+ }
+ else if (quadrature_name.compare("QGauss") == 0)
+ {
+ const std::pair<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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 <int dim,int spacedim>
+ FiniteElement<dim,spacedim> *get_fe_by_name (std::string &name)
+ {
+ return get_fe_by_name_ext<dim,spacedim> (name, fe_name_map[dim][spacedim]);
+ }
+ }
+ }
+
+
+
+ template <int dim>
+ FiniteElement<dim> *
+ get_fe_from_name (const std::string ¶meter_name)
+ {
+ return get_fe_by_name<dim,dim> (parameter_name);
+ }
+
+
+
+ template <int dim, int spacedim>
+ FiniteElement<dim, spacedim> *
+ 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<dim,spacedim> *fe = internal::get_fe_by_name<dim,spacedim> (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 <int dim, int spacedim>
+ void
+ compute_projection_from_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
+ const Quadrature<dim> &lhs_quadrature,
+ const Quadrature<dim> &rhs_quadrature,
+ FullMatrix<double> &X)
+ {
+ Assert (fe.n_components() == 1, ExcNotImplemented());
+
+ // first build the matrices M and Q
+ // described in the documentation
+ FullMatrix<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
+ FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
+
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
+ M(i,j) += fe.shape_value (i, lhs_quadrature.point(q)) *
+ fe.shape_value (j, lhs_quadrature.point(q)) *
+ lhs_quadrature.weight(q);
+
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
+ Q(i,q) += fe.shape_value (i, rhs_quadrature.point(q)) *
+ rhs_quadrature.weight(q);
+
+ // then invert M
+ FullMatrix<double> 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 <int dim, int spacedim>
+ void
+ compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
+ const Quadrature<dim> &quadrature,
+ FullMatrix<double> &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<quadrature.size(); ++q)
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ I_q(q,i) = fe.shape_value (i, quadrature.point(q));
+ }
+
+
+
+ template <int dim>
+ void
+ compute_projection_from_quadrature_points(
+ const FullMatrix<double> &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<double> component_at_node(n_support_points);
+ // component at the quadrature point
+ Vector<double> 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 <n_support_points; ++nn)
+ {
+ vector_of_tensors_at_nodes[nn][ii] = component_at_node(nn);
+ }
+ }
+ }
+
+
+
+ template <int dim>
+ void
+ compute_projection_from_quadrature_points(
+ const FullMatrix<double> &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<double> component_at_node(n_support_points);
+ // component at the quadrature point
+ Vector<double> 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 <n_support_points; ++nn)
+ {
+ (vector_of_tensors_at_nodes[nn])[row][column] = component_at_node(nn);
+ }
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ void compute_projection_from_face_quadrature_points_matrix(
+ const FiniteElement<dim, spacedim> &fe,
+ const Quadrature<dim - 1> &lhs_quadrature,
+ const Quadrature<dim - 1> &rhs_quadrature,
+ const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
+ const unsigned int face,
+ FullMatrix<double> &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<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
+ FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
+
+ {
+ // need an FEFaceValues object to evaluate shape function
+ // values on the specified face.
+ FEFaceValues <dim> 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.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
+ M(i,j) += fe_face_values.shape_value (i, q) *
+ fe_face_values.shape_value (j, q) *
+ lhs_quadrature.weight(q);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ M(i,i) = (M(i,i) == 0 ? 1 : M(i,i));
+ }
+ }
+
+ {
+ FEFaceValues <dim> 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<fe.dofs_per_cell; ++i)
+ for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
+ Q(i,q) += fe_face_values.shape_value (i, q) *
+ rhs_quadrature.weight(q);
+ }
+ // then invert M
+ FullMatrix<double> 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 <int dim>
+ void
+ hierarchic_to_lexicographic_numbering (unsigned int degree, std::vector<unsigned int> &h2l)
+ {
+ // number of support points in each direction
+ const unsigned int n = degree+1;
+
+ const unsigned int dofs_per_cell = Utilities::fixed_power<dim>(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<dofs_per_cell; ++i)
+ h2l[i] = i-1;
+
+ break;
+ }
+
+ case 2:
+ {
+ unsigned int next_index = 0;
+ // first the four vertices
+ h2l[next_index++] = 0;
+ h2l[next_index++] = n-1;
+ h2l[next_index++] = n*(n-1);
+ h2l[next_index++] = n*n-1;
+
+ // left line
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = (1+i)*n;
+
+ // right line
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = (2+i)*n-1;
+
+ // bottom line
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = 1+i;
+
+ // top line
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = n*(n-1)+i+1;
+
+ // inside quad
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ h2l[next_index++] = n*(i+1)+j+1;
+
+ Assert (next_index == dofs_per_cell, ExcInternalError());
+
+ break;
+ }
+
+ case 3:
+ {
+ unsigned int next_index = 0;
+ // first the eight vertices
+ h2l[next_index++] = 0; // 0
+ h2l[next_index++] = ( 1)*degree; // 1
+ h2l[next_index++] = ( n )*degree; // 2
+ h2l[next_index++] = ( n+1)*degree; // 3
+ h2l[next_index++] = (n*n )*degree; // 4
+ h2l[next_index++] = (n*n +1)*degree; // 5
+ h2l[next_index++] = (n*n+n )*degree; // 6
+ h2l[next_index++] = (n*n+n+1)*degree; // 7
+
+ // line 0
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = (i+1)*n;
+ // line 1
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = n-1+(i+1)*n;
+ // line 2
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = 1+i;
+ // line 3
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = 1+i+n*(n-1);
+
+ // line 4
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = (n-1)*n*n+(i+1)*n;
+ // line 5
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = (n-1)*(n*n+1)+(i+1)*n;
+ // line 6
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = n*n*(n-1)+i+1;
+ // line 7
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = n*n*(n-1)+i+1+n*(n-1);
+
+ // line 8
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = (i+1)*n*n;
+ // line 9
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = n-1+(i+1)*n*n;
+ // line 10
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = (i+1)*n*n+n*(n-1);
+ // line 11
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ h2l[next_index++] = n-1+(i+1)*n*n+n*(n-1);
+
+
+ // inside quads
+ // face 0
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ h2l[next_index++] = (i+1)*n*n+n*(j+1);
+ // face 1
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ h2l[next_index++] = (i+1)*n*n+n-1+n*(j+1);
+ // face 2, note the orientation!
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ h2l[next_index++] = (j+1)*n*n+i+1;
+ // face 3, note the orientation!
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ h2l[next_index++] = (j+1)*n*n+n*(n-1)+i+1;
+ // face 4
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ h2l[next_index++] = n*(i+1)+j+1;
+ // face 5
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ h2l[next_index++] = (n-1)*n*n+n*(i+1)+j+1;
+
+ // inside hex
+ for (unsigned int i=0; i<dofs_per_line; ++i)
+ for (unsigned int j=0; j<dofs_per_line; ++j)
+ for (unsigned int k=0; k<dofs_per_line; ++k)
+ h2l[next_index++] = n*n*(i+1)+n*(j+1)+k+1;
+
+ Assert (next_index == dofs_per_cell, ExcInternalError());
+
+ break;
+ }
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+
+
+
+ template <int dim>
+ void
+ hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &fe,
+ std::vector<unsigned int> &h2l)
+ {
+ Assert (h2l.size() == fe.dofs_per_cell,
+ ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell));
+ hierarchic_to_lexicographic_numbering<dim> (fe.dofs_per_line+1, h2l);
+ }
+
+
+
+ template <int dim>
+ std::vector<unsigned int>
+ hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &fe)
+ {
+ Assert (fe.n_components() == 1, ExcInvalidFE());
+ std::vector<unsigned int> h2l(fe.dofs_per_cell);
+ hierarchic_to_lexicographic_numbering<dim> (fe.dofs_per_line+1, h2l);
+ return (h2l);
+ }
+
+
+
+ template <int dim>
+ void
+ lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe,
+ std::vector<unsigned int> &l2h)
+ {
+ l2h = lexicographic_to_hierarchic_numbering (fe);
+ }
+
+
+
+ template <int dim>
+ std::vector<unsigned int>
+ lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe)
+ {
+ return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe));
+ }
+
+} // end of namespace FETools
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif /* dealii_fe_tools_templates_H */
// ---------------------------------------------------------------------
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/qprojector.h>
-#include <deal.II/base/thread_management.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/householder.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/fe/fe_face.h>
-#include <deal.II/fe/fe_q_iso_q1.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_q_dg0.h>
-#include <deal.II/fe/fe_q_bubbles.h>
-#include <deal.II/fe/fe_bernstein.h>
-#include <deal.II/fe/fe_q_hierarchical.h>
-#include <deal.II/fe/fe_dgq.h>
-#include <deal.II/fe/fe_dgp.h>
-#include <deal.II/fe/fe_dgp_monomial.h>
-#include <deal.II/fe/fe_dgp_nonparametric.h>
-#include <deal.II/fe/fe_dg_vector.h>
-#include <deal.II/fe/fe_nedelec.h>
-#include <deal.II/fe/fe_abf.h>
-#include <deal.II/fe/fe_bdm.h>
-#include <deal.II/fe/fe_raviart_thomas.h>
-#include <deal.II/fe/fe_rannacher_turek.h>
-#include <deal.II/fe/fe_nothing.h>
-#include <deal.II/fe/fe_system.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping_cartesian.h>
-#include <deal.II/fe/mapping_q1.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/hp/dof_handler.h>
-
-
-#include <deal.II/base/index_set.h>
-
-#include <cctype>
-#include <iostream>
-#include <memory>
-
+#include <deal.II/fe/fe_tools.templates.h>
DEAL_II_NAMESPACE_OPEN
-namespace FETools
-{
- namespace Compositing
- {
-
- template <int dim, int spacedim>
- FiniteElementData<dim>
- multiply_dof_numbers (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &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; i<fes.size(); i++)
- if (multiplicities[i]>0)
- {
- n_components = fes[i]->n_components();
- break;
- }
-
- for (unsigned int i=0; i<fes.size(); i++)
- if (multiplicities[i]>0)
- {
- 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<dim>::Conformity total_conformity
- = typename FiniteElementData<dim>::Conformity();
- {
- unsigned int index = 0;
- for (index=0; index<fes.size(); ++index)
- if (multiplicities[index]>0)
- {
- total_conformity = fes[index]->conforming_space;
- break;
- }
-
- for (; index<fes.size(); ++index)
- if (multiplicities[index]>0)
- total_conformity =
- typename FiniteElementData<dim>::Conformity(total_conformity
- &
- fes[index]->conforming_space);
- }
-
- std::vector<unsigned int> 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<dim> (dpo,
- (do_tensor_product ? multiplied_n_components : n_components),
- degree,
- total_conformity,
- block_indices);
- }
-
-
-
- template <int dim, int spacedim>
- FiniteElementData<dim>
- multiply_dof_numbers (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2,
- const unsigned int N2,
- const FiniteElement<dim,spacedim> *fe3,
- const unsigned int N3,
- const FiniteElement<dim,spacedim> *fe4,
- const unsigned int N4,
- const FiniteElement<dim,spacedim> *fe5,
- const unsigned int N5)
- {
- std::vector<const FiniteElement<dim,spacedim>*> fes;
- fes.push_back(fe1);
- fes.push_back(fe2);
- fes.push_back(fe3);
- fes.push_back(fe4);
- fes.push_back(fe5);
-
- std::vector<unsigned int> 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 <int dim, int spacedim>
- std::vector<bool>
- compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &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; i<fes.size(); ++i)
- if (multiplicities[i]>0) // 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<bool> 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<GeometryInfo<dim>::vertices_per_cell;
- ++vertex_number)
- {
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base]; ++m)
- 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 (index_in_base < fes[base]->dofs_per_cell,
- ExcInternalError());
- retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
- }
- }
-
- // 2. Lines
- if (GeometryInfo<dim>::lines_per_cell > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::lines_per_cell;
- ++line_number)
- {
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base]; ++m)
- 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 (index_in_base < fes[base]->dofs_per_cell,
- ExcInternalError());
- retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
- }
- }
-
- // 3. Quads
- if (GeometryInfo<dim>::quads_per_cell > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::quads_per_cell;
- ++quad_number)
- {
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base]; ++m)
- 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 (index_in_base < fes[base]->dofs_per_cell,
- ExcInternalError());
- retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
- }
- }
-
- // 4. Hexes
- if (GeometryInfo<dim>::hexes_per_cell > 0)
- for (unsigned int hex_number= 0;
- hex_number != GeometryInfo<dim>::hexes_per_cell;
- ++hex_number)
- {
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base]; ++m)
- 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 (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 <int dim, int spacedim>
- std::vector<bool>
- compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2,
- const unsigned int N2,
- const FiniteElement<dim,spacedim> *fe3,
- const unsigned int N3,
- const FiniteElement<dim,spacedim> *fe4,
- const unsigned int N4,
- const FiniteElement<dim,spacedim> *fe5,
- const unsigned int N5)
- {
- std::vector<const FiniteElement<dim,spacedim>*> fe_list;
- std::vector<unsigned int> 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 <int dim, int spacedim>
- std::vector<ComponentMask>
- compute_nonzero_components (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
- const std::vector<unsigned int> &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; i<fes.size(); ++i)
- if (multiplicities[i]>0) //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; i<fes.size(); ++i)
- if (multiplicities[i]>0) //needed because fe might be NULL
- n_components += fes[i]->n_components() * multiplicities[i];
- }
- else
- {
- for (unsigned int i=0; i<fes.size(); ++i)
- if (multiplicities[i]>0) //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; i<fes.size(); ++i)
- if (multiplicities[i]>0) //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<std::vector<bool> >
- retval (n_shape_functions, std::vector<bool> (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<GeometryInfo<dim>::vertices_per_cell;
- ++vertex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base];
- ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<dim>::lines_per_cell > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::lines_per_cell;
- ++line_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base];
- ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<dim>::quads_per_cell > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::quads_per_cell;
- ++quad_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base];
- ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<dim>::hexes_per_cell > 0)
- for (unsigned int hex_number= 0;
- hex_number != GeometryInfo<dim>::hexes_per_cell;
- ++hex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fes.size(); ++base)
- for (unsigned int m=0; m<multiplicities[base];
- ++m, comp_start+=fes[base]->n_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; c<fes[base]->n_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<vector<bool> > into a vector<ComponentMask>.
- // 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<bool> anyway.
- std::vector<ComponentMask> xretval (retval.size());
- for (unsigned int i=0; i<retval.size(); ++i)
- xretval[i] = ComponentMask(retval[i]);
- return xretval;
- }
-
-
-
- /**
- * Compute the non-zero vector components of a composed finite element.
- */
- template <int dim, int spacedim>
- std::vector<ComponentMask>
- compute_nonzero_components (const FiniteElement<dim,spacedim> *fe1,
- const unsigned int N1,
- const FiniteElement<dim,spacedim> *fe2,
- const unsigned int N2,
- const FiniteElement<dim,spacedim> *fe3,
- const unsigned int N3,
- const FiniteElement<dim,spacedim> *fe4,
- const unsigned int N4,
- const FiniteElement<dim,spacedim> *fe5,
- const unsigned int N5,
- const bool do_tensor_product)
- {
- std::vector<const FiniteElement<dim,spacedim>*> fe_list;
- std::vector<unsigned int> 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 <int dim, int spacedim>
- 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<dim,spacedim> &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<fe.base_element(base).n_components(); ++k)
- component_to_base_table[total_index++]
- = std::make_pair(std::make_pair(base,k), m);
- }
- Assert (total_index == component_to_base_table.size(),
- ExcInternalError());
- }
- else
- {
- // The base element establishing a component does not make sense in this case.
- // Set up to something meaningless:
- for (unsigned int i = 0; i < component_to_base_table.size(); i++)
- component_to_base_table[i] = std::make_pair(std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int), numbers::invalid_unsigned_int);
-
- }
-
-
- // Initialize index tables. Multi-component base elements have to be
- // thought of. For non-primitive shape functions, have a special invalid
- // index.
- const std::pair<unsigned int, unsigned int>
- 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<GeometryInfo<dim>::vertices_per_cell;
- ++vertex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_vertex;
- ++local_index, ++total_index)
- {
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_vertex*vertex_number +
- local_index);
-
- system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base, m), index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).system_to_component_index(index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int index_in_comp
- = fe.base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
- = std::make_pair (comp, index_in_comp);
- }
- else
- system_to_component_table[total_index] = non_primitive_index;
- }
- }
-
- // 2. Lines
- if (GeometryInfo<dim>::lines_per_cell > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::lines_per_cell;
- ++line_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_line;
- ++local_index, ++total_index)
- {
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_line*line_number +
- local_index +
- fe.base_element(base).first_line_index);
-
- system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).system_to_component_index(index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int index_in_comp
- = fe.base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
- = std::make_pair (comp, index_in_comp);
- }
- else
- system_to_component_table[total_index] = non_primitive_index;
- }
- }
-
- // 3. Quads
- if (GeometryInfo<dim>::quads_per_cell > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::quads_per_cell;
- ++quad_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_quad;
- ++local_index, ++total_index)
- {
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_quad*quad_number +
- local_index +
- fe.base_element(base).first_quad_index);
-
- system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).system_to_component_index(index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int index_in_comp
- = fe.base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
- = std::make_pair (comp, index_in_comp);
- }
- else
- system_to_component_table[total_index] = non_primitive_index;
- }
- }
-
- // 4. Hexes
- if (GeometryInfo<dim>::hexes_per_cell > 0)
- for (unsigned int hex_number= 0;
- hex_number != GeometryInfo<dim>::hexes_per_cell;
- ++hex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start+=fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_hex;
- ++local_index, ++total_index)
- {
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_hex*hex_number +
- local_index +
- fe.base_element(base).first_hex_index);
-
- system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).system_to_component_index(index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int index_in_comp
- = fe.base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
- = std::make_pair (comp, index_in_comp);
- }
- else
- system_to_component_table[total_index] = non_primitive_index;
- }
- }
- }
-
-
-
- template <int dim, int spacedim>
- 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<dim,spacedim> &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<unsigned int, unsigned int>
- 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<GeometryInfo<dim>::vertices_per_face;
- ++vertex_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_vertex;
- ++local_index, ++total_index)
- {
- // get (cell) index of this shape function inside the base
- // element to see whether the shape function is primitive
- // (assume that all shape functions on vertices share the same
- // primitivity property; assume likewise for all shape functions
- // located on lines, quads, etc. this way, we can ask for
- // primitivity of only _one_ shape function, which is taken as
- // representative for all others located on the same type of
- // object):
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_vertex*vertex_number +
- local_index);
-
- const unsigned int face_index_in_base
- = (fe.base_element(base).dofs_per_vertex*vertex_number +
- local_index);
-
- face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int face_index_in_comp
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
- }
- else
- face_system_to_component_table[total_index] = non_primitive_index;
- }
- }
-
- // 2. Lines
- if (GeometryInfo<dim>::lines_per_face > 0)
- for (unsigned int line_number= 0;
- line_number != GeometryInfo<dim>::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<fe.element_multiplicity(base);
- ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_line;
- ++local_index, ++total_index)
- {
- // do everything alike for this type of object
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_line*line_number +
- local_index +
- fe.base_element(base).first_line_index);
-
- const unsigned int face_index_in_base
- = (fe.base_element(base).first_face_line_index +
- fe.base_element(base).dofs_per_line * line_number +
- local_index);
-
- face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int face_index_in_comp
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
- }
- else
- face_system_to_component_table[total_index] = non_primitive_index;
- }
- }
-
- // 3. Quads
- if (GeometryInfo<dim>::quads_per_face > 0)
- for (unsigned int quad_number= 0;
- quad_number != GeometryInfo<dim>::quads_per_face;
- ++quad_number)
- {
- unsigned int comp_start = 0;
- for (unsigned int base=0; base<fe.n_base_elements(); ++base)
- for (unsigned int m=0; m<fe.element_multiplicity(base);
- ++m, comp_start += fe.base_element(base).n_components() * do_tensor_product)
- for (unsigned int local_index = 0;
- local_index < fe.base_element(base).dofs_per_quad;
- ++local_index, ++total_index)
- {
- // do everything alike for this type of object
- const unsigned int index_in_base
- = (fe.base_element(base).dofs_per_quad*quad_number +
- local_index +
- fe.base_element(base).first_quad_index);
-
- const unsigned int face_index_in_base
- = (fe.base_element(base).first_face_quad_index +
- fe.base_element(base).dofs_per_quad * quad_number +
- local_index);
-
- face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
- if (fe.base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).first;
- const unsigned int comp
- = comp_start + comp_in_base;
- const unsigned int face_index_in_comp
- = fe.base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
- }
- else
- face_system_to_component_table[total_index] = non_primitive_index;
- }
- }
- Assert (total_index == fe.dofs_per_face, ExcInternalError());
- Assert (total_index == face_system_to_component_table.size(),
- ExcInternalError());
- Assert (total_index == face_system_to_base_table.size(),
- ExcInternalError());
- }
- }
-
-
-
- // Not implemented in the general case.
- template <class FE>
- FiniteElement<FE::dimension, FE::space_dimension> *
- FEFactory<FE>::get (const Quadrature<1> &) const
- {
- Assert(false, ExcNotImplemented());
- return nullptr;
- }
-
- // Specializations for FE_Q.
- template <>
- FiniteElement<1, 1> *
- FEFactory<FE_Q<1, 1> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q<1>(quad);
- }
-
- template <>
- FiniteElement<2, 2> *
- FEFactory<FE_Q<2, 2> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q<2>(quad);
- }
-
- template <>
- FiniteElement<3, 3> *
- FEFactory<FE_Q<3, 3> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q<3>(quad);
- }
-
- // Specializations for FE_Q_DG0.
- template <>
- FiniteElement<1, 1> *
- FEFactory<FE_Q_DG0<1, 1> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q_DG0<1>(quad);
- }
-
- template <>
- FiniteElement<2, 2> *
- FEFactory<FE_Q_DG0<2, 2> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q_DG0<2>(quad);
- }
-
- template <>
- FiniteElement<3, 3> *
- FEFactory<FE_Q_DG0<3, 3> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q_DG0<3>(quad);
- }
-
- // Specializations for FE_Q_Bubbles.
- template <>
- FiniteElement<1, 1> *
- FEFactory<FE_Q_Bubbles<1, 1> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q_Bubbles<1>(quad);
- }
-
- template <>
- FiniteElement<2, 2> *
- FEFactory<FE_Q_Bubbles<2, 2> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q_Bubbles<2>(quad);
- }
-
- template <>
- FiniteElement<3, 3> *
- FEFactory<FE_Q_Bubbles<3, 3> >::get (const Quadrature<1> &quad) const
- {
- return new FE_Q_Bubbles<3>(quad);
- }
-
- // Specializations for FE_DGQArbitraryNodes.
- template <>
- FiniteElement<1, 1> *
- FEFactory<FE_DGQ<1> >::get (const Quadrature<1> &quad) const
- {
- return new FE_DGQArbitraryNodes<1>(quad);
- }
-
- template <>
- FiniteElement<1, 2> *
- FEFactory<FE_DGQ<1, 2> >::get (const Quadrature<1> &quad) const
- {
- return new FE_DGQArbitraryNodes<1, 2>(quad);
- }
-
- template <>
- FiniteElement<1, 3> *
- FEFactory<FE_DGQ<1, 3> >::get (const Quadrature<1> &quad) const
- {
- return new FE_DGQArbitraryNodes<1, 3>(quad);
- }
-
- template <>
- FiniteElement<2, 2> *
- FEFactory<FE_DGQ<2> >::get (const Quadrature<1> &quad) const
- {
- return new FE_DGQArbitraryNodes<2>(quad);
- }
-
- template <>
- FiniteElement<2, 3> *
- FEFactory<FE_DGQ<2, 3> >::get (const Quadrature<1> &quad) const
- {
- return new FE_DGQArbitraryNodes<2, 3>(quad);
- }
-
- template <>
- FiniteElement<3, 3> *
- FEFactory<FE_DGQ<3> >::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 <int dim>
- void
- fill_no_codim_fe_names (std::map<std::string,std::shared_ptr<const Subscriptor> > &result)
- {
- typedef std::shared_ptr<const Subscriptor> FEFactoryPointer;
-
- result["FE_Q_Hierarchical"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q_Hierarchical<dim> >);
- result["FE_ABF"]
- = FEFactoryPointer(new FETools::FEFactory<FE_ABF<dim> >);
- result["FE_Bernstein"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Bernstein<dim> >);
- result["FE_BDM"]
- = FEFactoryPointer(new FETools::FEFactory<FE_BDM<dim> >);
- result["FE_DGBDM"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGBDM<dim> >);
- result["FE_DGNedelec"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGNedelec<dim> >);
- result["FE_DGRaviartThomas"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGRaviartThomas<dim> >);
- result["FE_RaviartThomas"]
- = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomas<dim> >);
- result["FE_RaviartThomasNodal"]
- = FEFactoryPointer(new FETools::FEFactory<FE_RaviartThomasNodal<dim> >);
- result["FE_Nedelec"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Nedelec<dim> >);
- result["FE_DGPNonparametric"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGPNonparametric<dim> >);
- result["FE_DGP"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGP<dim> >);
- result["FE_DGPMonomial"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGPMonomial<dim> >);
- result["FE_DGQ"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim> >);
- result["FE_DGQArbitraryNodes"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim> >);
- result["FE_DGQLegendre"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQLegendre<dim> >);
- result["FE_DGQHermite"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQHermite<dim> >);
- result["FE_FaceQ"]
- = FEFactoryPointer(new FETools::FEFactory<FE_FaceQ<dim> >);
- result["FE_FaceP"]
- = FEFactoryPointer(new FETools::FEFactory<FE_FaceP<dim> >);
- result["FE_Q"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q<dim> >);
- result["FE_Q_DG0"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q_DG0<dim> >);
- result["FE_Q_Bubbles"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q_Bubbles<dim> >);
- result["FE_Q_iso_Q1"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q_iso_Q1<dim> >);
- result["FE_Nothing"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Nothing<dim> >);
- result["FE_RannacherTurek"]
- = FEFactoryPointer(new FETools::FEFactory<FE_RannacherTurek<dim> >);
- }
-
-
-
- // This function fills a map from names to finite elements for any
- // dimension and codimension for those elements which support
- // nonzero codimension.
- template <int dim, int spacedim>
- void
- fill_codim_fe_names (std::map<std::string,std::shared_ptr<const Subscriptor> > &result)
- {
- typedef std::shared_ptr<const Subscriptor> FEFactoryPointer;
-
- result["FE_Bernstein"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Bernstein<dim,spacedim> >);
- result["FE_DGP"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGP<dim,spacedim> >);
- result["FE_DGQ"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim,spacedim> >);
- result["FE_Nothing"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Nothing<dim,spacedim> >);
- result["FE_DGQArbitraryNodes"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQ<dim,spacedim> >);
- result["FE_DGQLegendre"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQLegendre<dim,spacedim> >);
- result["FE_DGQHermite"]
- = FEFactoryPointer(new FETools::FEFactory<FE_DGQHermite<dim,spacedim> >);
- result["FE_Q_Bubbles"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q_Bubbles<dim,spacedim> >);
- result["FE_Q_DG0"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q_DG0<dim,spacedim> >);
- result["FE_Q_iso_Q1"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q_iso_Q1<dim,spacedim> >);
- result["FE_Q"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Q<dim,spacedim> >);
- result["FE_Bernstein"]
- = FEFactoryPointer(new FETools::FEFactory<FE_Bernstein<dim,spacedim> >);
- }
-
- // 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<std::vector<
- std::map<std::string,
- std::shared_ptr<const Subscriptor> > > >
- fill_default_map()
- {
- std::vector<std::vector<
- std::map<std::string,
- std::shared_ptr<const Subscriptor> > > >
- 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<dim,spacedim> 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<std::vector<
- std::map<std::string,
- std::shared_ptr<const Subscriptor> > > >
- fe_name_map = fill_default_map();
-}
-
-
-
-
-
-
-namespace
-{
-
- // forwarder function for
- // FE::get_interpolation_matrix. we
- // will want to call that function
- // for arbitrary FullMatrix<T>
- // 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 <int dim, int spacedim>
- inline
- void gim_forwarder (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<double> &interpolation_matrix)
- {
- fe2.get_interpolation_matrix (fe1, interpolation_matrix);
- }
-
-
-
- template <int dim, typename number, int spacedim>
- inline
- void gim_forwarder (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &interpolation_matrix)
- {
- FullMatrix<double> 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 "<dim>" or the
- // specialized string with "dim"
- // replaced with the numeric value
- // of the template argument
- template <int dim, int spacedim>
- 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 <int dim, int spacedim>
- FEFactoryBase<dim,spacedim>::~FEFactoryBase()
- {}
-
-
-
- template <int dim, int spacedim>
- void compute_component_wise(
- const FiniteElement<dim,spacedim> &element,
- std::vector<unsigned int> &renumbering,
- std::vector<std::vector<unsigned int> > &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<comp_start.size(); ++i)
- {
- comp_start[i].resize(element.element_multiplicity(i));
- const unsigned int increment
- = element.base_element(i).dofs_per_cell;
-
- for (unsigned int j=0; j<comp_start[i].size(); ++j)
- {
- comp_start[i][j] = k;
- k += increment;
- }
- }
-
- // For each index i of the
- // unstructured cellwise
- // numbering, renumbering
- // contains the index of the
- // cell-block numbering
- for (unsigned int i=0; i<element.dofs_per_cell; ++i)
- {
- std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
- indices = element.system_to_base_index(i);
- renumbering[i] = comp_start[indices.first.first][indices.first.second]
- +indices.second;
- }
- }
-
-
-
- template <int dim, int spacedim>
- void compute_block_renumbering (
- const FiniteElement<dim,spacedim> &element,
- std::vector<types::global_dof_index> &renumbering,
- std::vector<types::global_dof_index> &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<element.n_base_elements(); ++b)
- for (unsigned int m=0; m<element.element_multiplicity(b); ++m)
- {
- block_data[count++] = (return_start_indices)
- ? k
- : (element.base_element(b).n_dofs_per_cell());
- k += element.base_element(b).n_dofs_per_cell();
- }
- Assert (count == element.n_blocks(), ExcInternalError());
-
- std::vector<types::global_dof_index> start_indices(block_data.size());
- k = 0;
- for (unsigned int i=0; i<block_data.size(); ++i)
- if (return_start_indices)
- start_indices[i] = block_data[i];
- else
- {
- start_indices[i] = k;
- k += block_data[i];
- }
-
- for (unsigned int i=0; i<element.dofs_per_cell; ++i)
- {
- std::pair<unsigned int, types::global_dof_index>
- indices = element.system_to_block_index(i);
- renumbering[i] = start_indices[indices.first]
- +indices.second;
- }
- }
-
-
-
- template <int dim, typename number, int spacedim>
- void get_interpolation_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<dim,spacedim>::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<Point<dim> > &
- fe2_support_points = fe2.get_unit_support_points ();
-
- typedef FiniteElement<dim,spacedim> FEL;
- Assert(fe2_support_points.size()==fe2.dofs_per_cell,
- typename FEL::ExcFEHasNoSupportPoints());
-
- for (unsigned int i=0; i<fe2.dofs_per_cell; ++i)
- {
- const unsigned int i1 = fe2.system_to_component_index(i).first;
- for (unsigned int j=0; j<fe1.dofs_per_cell; ++j)
- {
- const unsigned int j1 = fe1.system_to_component_index(j).first;
- if (i1==j1)
- interpolation_matrix(i,j) = fe1.shape_value (j,fe2_support_points[i]);
- else
- interpolation_matrix(i,j)=0.;
- }
- }
- }
-
-
-
- template <int dim, typename number, int spacedim>
- void get_back_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<number> first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell);
- FullMatrix<number> 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 <int dim, typename number, int spacedim>
- void get_interpolation_difference_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<number> interpolation_matrix(fe1.dofs_per_cell);
- get_back_interpolation_matrix(fe1, fe2, interpolation_matrix);
-
- for (unsigned int i=0; i<fe1.dofs_per_cell; ++i)
- difference_matrix(i,i) = 1.;
-
- // compute difference
- difference_matrix.add (-1, interpolation_matrix);
- }
-
-
-
- template <int dim, typename number, int spacedim>
- void get_projection_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<dim,spacedim> 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<dim> quadrature(degree+1);
- // Set up FEValues.
- const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values;
- FEValues<dim> val1 (fe1, quadrature, update_values);
- val1.reinit (tr.begin_active());
- FEValues<dim> val2 (fe2, quadrature, flags);
- val2.reinit (tr.begin_active());
-
- // Integrate and invert mass matrix
- // This happens in the target space
- FullMatrix<double> mass (n2, n2);
-
- for (unsigned int k=0; k<quadrature.size(); ++k)
- {
- const double w = val2.JxW(k);
- for (unsigned int i=0; i<n2; ++i)
- {
- const double v = val2.shape_value(i,k);
- for (unsigned int j=0; j<n2; ++j)
- mass(i,j) += w*v * val2.shape_value(j,k);
- }
- }
- // Gauss-Jordan should be
- // sufficient since we expect the
- // mass matrix to be
- // well-conditioned
- mass.gauss_jordan();
-
- // Now, test every function of fe1
- // with test functions of fe2 and
- // compute the projection of each
- // unit vector.
- Vector<double> b(n2);
- Vector<double> x(n2);
-
- for (unsigned int j=0; j<n1; ++j)
- {
- b = 0.;
- for (unsigned int i=0; i<n2; ++i)
- for (unsigned int k=0; k<quadrature.size(); ++k)
- {
- const double w = val2.JxW(k);
- const double u = val1.shape_value(j,k);
- const double v = val2.shape_value(i,k);
- b(i) += u*v*w;
- }
-
- // Multiply by the inverse
- mass.vmult(x,b);
- for (unsigned int i=0; i<n2; ++i)
- matrix(i,j) = x(i);
- }
- }
-
-
-
- template <int dim, int spacedim>
- FullMatrix<double>
- compute_node_matrix(const FiniteElement<dim,spacedim> &fe)
- {
- const unsigned int n_dofs = fe.dofs_per_cell;
-
- FullMatrix<double> N (n_dofs, n_dofs);
-
- Assert (fe.has_generalized_support_points(), ExcNotInitialized());
- Assert (fe.n_components() == dim, ExcNotImplemented());
-
- const std::vector<Point<dim> > &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<Vector<double> >
- support_point_values (points.size(), Vector<double>(dim));
-
- // In this vector, we store the
- // result of the interpolation
- std::vector<double> 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<n_dofs; ++i)
- {
- // get the values of the current set of shape functions
- // at the generalized support points
- for (unsigned int k=0; k<points.size(); ++k)
- for (unsigned int d=0; d<dim; ++d)
- {
- support_point_values[k][d] = fe.shape_value_component(i, points[k], d);
- Assert (numbers::is_finite(support_point_values[k][d]), ExcInternalError());
- }
-
- fe.convert_generalized_support_point_values_to_dof_values(support_point_values,
- nodal_values);
-
- // Enter the interpolated dofs into the matrix
- for (unsigned int j=0; j<n_dofs; ++j)
- {
- N(j,i) = nodal_values[j];
- Assert (numbers::is_finite(nodal_values[j]), ExcInternalError());
- }
- }
-
- return N;
- }
-
-
-
- /*
- template <>
- void
- compute_embedding_matrices(const FiniteElement<1,2> &,
- std::vector<std::vector<FullMatrix<double> > > &,
- const bool)
- {
- Assert(false, ExcNotImplemented());
- }
-
-
- template <>
- void
- compute_embedding_matrices(const FiniteElement<1,3> &,
- std::vector<std::vector<FullMatrix<double> > > &,
- const bool)
- {
- Assert(false, ExcNotImplemented());
- }
-
-
-
- template <>
- void
- compute_embedding_matrices(const FiniteElement<2,3>&,
- std::vector<std::vector<FullMatrix<double> > >&,
- const bool)
- {
- Assert(false, ExcNotImplemented());
- }
-
- */
-
- namespace
- {
- template <int dim, typename number, int spacedim>
- void
- compute_embedding_for_shape_function (
- const unsigned int i,
- const FiniteElement<dim, spacedim> &fe,
- const FEValues<dim, spacedim> &coarse,
- const Householder<double> &H,
- FullMatrix<number> &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<number> v_coarse(nq*nd);
- Vector<number> 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 <int dim, typename number, int spacedim>
- void
- compute_embedding_matrices_for_refinement_case (
- const FiniteElement<dim, spacedim> &fe,
- std::vector<FullMatrix<number> > &matrices,
- const unsigned int ref_case,
- const double threshold)
- {
- const unsigned int n = fe.dofs_per_cell;
- const unsigned int nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(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<dim,spacedim> tria;
- GridGenerator::hyper_cube (tria, 0, 1);
- tria.begin_active()->set_refine_flag (RefinementCase<dim>(ref_case));
- tria.execute_coarsening_and_refinement ();
-
- const unsigned int degree = fe.degree;
- QGauss<dim> q_fine (degree+1);
- const unsigned int nq = q_fine.size();
-
- FEValues<dim,spacedim> 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<number> 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<double> H (A);
- unsigned int cell_number = 0;
-
- Threads::TaskGroup<void> task_group;
-
- for (typename Triangulation<dim,spacedim>::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<Point<spacedim> > &q_points_fine = fine.get_quadrature_points();
- std::vector<Point<dim> > q_points_coarse(q_points_fine.size());
- for (unsigned int i=0; i<q_points_fine.size(); ++i)
- for (unsigned int j=0; j<dim; ++j)
- q_points_coarse[i](j) = q_points_fine[i](j);
- const Quadrature<dim> q_coarse (q_points_coarse,
- fine.get_JxW_values ());
- FEValues<dim,spacedim> coarse (fe, q_coarse, update_values);
-
- coarse.reinit (tria.begin (0));
-
- FullMatrix<double> &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<dim, number, spacedim>,
- 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<dim, number, spacedim>
- (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<dim>::n_children (RefinementCase<dim> (ref_case)),
- ExcInternalError ());
- }
- }
-
-
-
- template <int dim, typename number, int spacedim>
- void
- compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
- std::vector<std::vector<FullMatrix<number> > > &matrices,
- const bool isotropic_only,
- const double threshold)
- {
- Threads::TaskGroup<void> task_group;
-
- // loop over all possible refinement cases
- unsigned int ref_case = (isotropic_only)
- ? RefinementCase<dim>::isotropic_refinement
- : RefinementCase<dim>::cut_x;
-
- for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
- task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case<dim, number, spacedim>,
- fe, matrices[ref_case-1], ref_case, threshold);
-
- task_group.join_all ();
- }
-
-
-
- template <int dim, typename number, int spacedim>
- void
- compute_face_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
- FullMatrix<number> (&matrices)[GeometryInfo<dim>::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<dim>::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<dim>::Hdiv);
- const bool tangential = fe.conforms(FiniteElementData<dim>::Hcurl);
-
- 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));
- }
-
- // In order to make the loops below
- // simpler, we introduce vectors
- // containing for indices 0-n the
- // number of the corresponding
- // shape value on the cell.
- std::vector<unsigned int> face_c_dofs(n);
- std::vector<unsigned int> face_f_dofs(n);
- {
- unsigned int face_dof=0;
- for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
- {
- const unsigned int offset_c = GeometryInfo<dim>::face_to_cell_vertices(face_coarse, i)
- *fe.dofs_per_vertex;
- const unsigned int offset_f = GeometryInfo<dim>::face_to_cell_vertices(face_fine, i)
- *fe.dofs_per_vertex;
- for (unsigned int j=0; j<fe.dofs_per_vertex; ++j)
- {
- face_c_dofs[face_dof] = offset_c + j;
- face_f_dofs[face_dof] = offset_f + j;
- ++face_dof;
- }
- }
- for (unsigned int i=1; i<=GeometryInfo<dim>::lines_per_face; ++i)
- {
- const unsigned int offset_c = fe.first_line_index
- + GeometryInfo<dim>::face_to_cell_lines(face_coarse, i-1)
- *fe.dofs_per_line;
- const unsigned int offset_f = fe.first_line_index
- + GeometryInfo<dim>::face_to_cell_lines(face_fine, i-1)
- *fe.dofs_per_line;
- for (unsigned int j=0; j<fe.dofs_per_line; ++j)
- {
- face_c_dofs[face_dof] = offset_c + j;
- face_f_dofs[face_dof] = offset_f + j;
- ++face_dof;
- }
- }
- for (unsigned int i=1; i<=GeometryInfo<dim>::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<fe.dofs_per_quad; ++j)
- {
- face_c_dofs[face_dof] = offset_c + j;
- face_f_dofs[face_dof] = offset_f + j;
- ++face_dof;
- }
- }
- Assert (face_dof == fe.dofs_per_face, ExcInternalError());
- }
-
- // Set up meshes, one with a single
- // reference cell and refine it once
- Triangulation<dim,spacedim> tria;
- GridGenerator::hyper_cube (tria, 0, 1);
- tria.refine_global(1);
- MappingCartesian<dim> 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<dim-1> q_gauss(degree+1);
- const Quadrature<dim> q_fine = QProjector<dim>::project_to_face(q_gauss, face_fine);
- const unsigned int nq = q_fine.size();
-
- FEValues<dim> 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<number> A(nq*nd, n);
- for (unsigned int j=0; j<n; ++j)
- for (unsigned int k=0; k<nq; ++k)
- if (nd != dim)
- for (unsigned int d=0; d<nd; ++d)
- A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
- else
- {
- if (normal)
- A(k*nd,j) = fine.shape_value_component(face_f_dofs[j],k,0);
- if (tangential)
- for (unsigned int d=1; d<dim; ++d)
- A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
- }
-
- Householder<double> H(A);
-
- Vector<number> v_coarse(nq*nd);
- Vector<number> v_fine(n);
-
-
-
- for (unsigned int cell_number = 0; cell_number < GeometryInfo<dim>::max_children_per_face;
- ++cell_number)
- {
- const Quadrature<dim> q_coarse
- = QProjector<dim>::project_to_subface(q_gauss, face_coarse, cell_number);
- FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
-
- typename Triangulation<dim,spacedim>::active_cell_iterator fine_cell
- = tria.begin(0)->child(GeometryInfo<dim>::child_cell_on_face(
- tria.begin(0)->refinement_case(), face_coarse, cell_number));
- fine.reinit(fine_cell);
- coarse.reinit(tria.begin(0));
-
- FullMatrix<double> &this_matrix = matrices[cell_number];
-
- // Compute this once for each
- // coarse grid basis function
- for (unsigned int i=0; i<n; ++i)
- {
- // The right hand side of
- // the least squares
- // problem consists of the
- // function values of the
- // coarse grid function in
- // each quadrature point.
- for (unsigned int k=0; k<nq; ++k)
- if (nd != dim)
- for (unsigned int d=0; d<nd; ++d)
- v_coarse(k*nd+d) = coarse.shape_value_component (face_c_dofs[i],k,d);
- else
- {
- if (normal)
- v_coarse(k*nd) = coarse.shape_value_component(face_c_dofs[i],k,0);
- if (tangential)
- for (unsigned int d=1; d<dim; ++d)
- v_coarse(k*nd+d) = coarse.shape_value_component(face_c_dofs[i],k,d);
- }
- // solve the least squares
- // problem.
- const double result = H.least_squares(v_fine, v_coarse);
- Assert (result <= threshold, ExcLeastSquaresError(result));
- // Avoid compiler 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);
- }
- // 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.;
- }
- }
-
-
-
- template <int dim, typename number, int spacedim>
- void
- compute_projection_matrices(const FiniteElement<dim,spacedim> &fe,
- std::vector<std::vector<FullMatrix<number> > > &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<dim> q_fine(degree+1);
- const unsigned int nq = q_fine.size();
-
- // create mass matrix on coarse cell.
- FullMatrix<number> mass(n, n);
- {
- // set up a triangulation for coarse cell
- Triangulation<dim,spacedim> tr;
- GridGenerator::hyper_cube (tr, 0, 1);
-
- FEValues<dim,spacedim> coarse (fe, q_fine,
- update_JxW_values | update_values);
-
- typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
- = tr.begin(0);
- coarse.reinit (coarse_cell);
-
- const std::vector<double> &JxW = coarse.get_JxW_values();
- for (unsigned int i=0; i<n; ++i)
- for (unsigned int j=0; j<n; ++j)
- if (fe.is_primitive())
- {
- const double *coarse_i = &coarse.shape_value(i,0);
- const double *coarse_j = &coarse.shape_value(j,0);
- double mass_ij = 0;
- for (unsigned int k=0; k<nq; ++k)
- mass_ij += JxW[k] * coarse_i[k] * coarse_j[k];
- mass(i,j) = mass_ij;
- }
- else
- {
- double mass_ij = 0;
- for (unsigned int d=0; d<nd; ++d)
- for (unsigned int k=0; k<nq; ++k)
- mass_ij += JxW[k] * coarse.shape_value_component(i,k,d)
- * coarse.shape_value_component(j,k,d);
- mass(i,j) = mass_ij;
- }
-
- // invert mass matrix
- mass.gauss_jordan();
- }
-
- // loop over all possible
- // refinement cases
- unsigned int ref_case = (isotropic_only)
- ? RefinementCase<dim>::isotropic_refinement
- : RefinementCase<dim>::cut_x;
- for (; ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
- {
- const unsigned int
- nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
-
- for (unsigned int i=0; i<nc; ++i)
- {
- Assert(matrices[ref_case-1][i].n() == n,
- ExcDimensionMismatch(matrices[ref_case-1][i].n(),n));
- Assert(matrices[ref_case-1][i].m() == n,
- ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
- }
-
- // create a respective refinement on the
- // triangulation
- Triangulation<dim,spacedim> tr;
- GridGenerator::hyper_cube (tr, 0, 1);
- tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
- tr.execute_coarsening_and_refinement();
-
- FEValues<dim,spacedim> fine (StaticMappingQ1<dim,spacedim>::mapping, fe, q_fine,
- update_quadrature_points | update_JxW_values |
- update_values);
-
- typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
- = tr.begin(0);
-
- Vector<number> v_coarse(n);
- Vector<number> v_fine(n);
-
- for (unsigned int cell_number=0; cell_number<nc; ++cell_number)
- {
- FullMatrix<double> &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<Point<spacedim> > &q_points_fine = fine.get_quadrature_points();
- std::vector<Point<dim> > q_points_coarse(q_points_fine.size());
- for (unsigned int q=0; q<q_points_fine.size(); ++q)
- for (unsigned int j=0; j<dim; ++j)
- q_points_coarse[q](j) = q_points_fine[q](j);
- Quadrature<dim> q_coarse (q_points_coarse,
- fine.get_JxW_values());
- FEValues<dim,spacedim> coarse (StaticMappingQ1<dim,spacedim>::mapping, fe, q_coarse, update_values);
- coarse.reinit(coarse_cell);
-
- // Build RHS
-
- const std::vector<double> &JxW = fine.get_JxW_values();
-
- // Outer loop over all fine
- // grid shape functions phi_j
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- {
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- {
- if (fe.is_primitive())
- {
- const double *coarse_i = &coarse.shape_value(i,0);
- const double *fine_j = &fine.shape_value(j,0);
-
- double update = 0;
- for (unsigned int k=0; k<nq; ++k)
- update += JxW[k] * coarse_i[k] * fine_j[k];
- v_fine(i) = update;
- }
- else
- {
- double update = 0;
- for (unsigned int d=0; d<nd; ++d)
- for (unsigned int k=0; k<nq; ++k)
- update += JxW[k] * coarse.shape_value_component(i,k,d)
- * fine.shape_value_component(j,k,d);
- v_fine(i) = update;
- }
- }
-
- // RHS ready. Solve system
- // and enter row into
- // matrix
- mass.vmult (v_coarse, v_fine);
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- this_matrix(i,j) = v_coarse(i);
- }
-
- // 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.;
- }
- }
- }
-
-
-
- template <int dim, int spacedim>
- void
- add_fe_name(const std::string ¶meter_name,
- const FEFactoryBase<dim,spacedim> *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<const Subscriptor> (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 <int dim, int spacedim>
- FiniteElement<dim,spacedim> *
- get_fe_by_name_ext (std::string &name,
- const std::map<std::string,
- std::shared_ptr<const Subscriptor> >
- &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<const FiniteElement<dim,spacedim>*> base_fes;
- std::vector<unsigned int> 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<dim,spacedim> (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<int,unsigned int> 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<dim,spacedim> *system_element = nullptr;
-
- // uses new FESystem constructor
- // which is independent of
- // the number of FEs in the system
- system_element = new FESystem<dim,spacedim>(base_fes, base_multiplicities);
-
- // now we don't need the
- // list of base elements
- // any more
- for (unsigned int i=0; i<base_fes.size(); ++i)
- delete base_fes[i];
-
- // finally return our
- // findings
- // Add the closing ']' to
- // the length
- return system_element;
-
- }
- catch (...)
- {
- // ups, some exception
- // was thrown. prevent a
- // memory leak, and then
- // pass on the exception
- // to the caller
- for (unsigned int i=0; i<base_fes.size(); ++i)
- delete base_fes[i];
- throw;
- }
-
- // this is a place where we
- // should really never get,
- // since above we have either
- // returned from the
- // try-clause, or have
- // re-thrown in the catch
- // clause. check that we
- // never get here
- Assert (false, ExcInternalError());
- }
- else if (name_part == "FE_Nothing")
- {
- // remove the () from FE_Nothing()
- name.erase(0,2);
-
- // this is a bit of a hack, as
- // FE_Nothing does not take a
- // degree, but it does take an
- // argument, which defaults to 1,
- // so this properly returns
- // FE_Nothing()
- const Subscriptor *ptr = fe_name_map.find(name_part)->second.get();
- const FEFactoryBase<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(ptr);
- return fef->get(QGaussLobatto<1>(tmp.first));
- }
- else if (quadrature_name.compare("QGauss") == 0)
- {
- const std::pair<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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<int,unsigned int> 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<dim,spacedim> *fef=dynamic_cast<const FEFactoryBase<dim,spacedim>*>(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 <int dim,int spacedim>
- FiniteElement<dim,spacedim> *get_fe_by_name (std::string &name)
- {
- return get_fe_by_name_ext<dim,spacedim> (name, fe_name_map[dim][spacedim]);
- }
- }
- }
-
-
-
- template <int dim>
- FiniteElement<dim> *
- get_fe_from_name (const std::string ¶meter_name)
- {
- return get_fe_by_name<dim,dim> (parameter_name);
- }
-
-
-
- template <int dim, int spacedim>
- FiniteElement<dim, spacedim> *
- 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<dim,spacedim> *fe = internal::get_fe_by_name<dim,spacedim> (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 <int dim, int spacedim>
- void
- compute_projection_from_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
- const Quadrature<dim> &lhs_quadrature,
- const Quadrature<dim> &rhs_quadrature,
- FullMatrix<double> &X)
- {
- Assert (fe.n_components() == 1, ExcNotImplemented());
-
- // first build the matrices M and Q
- // described in the documentation
- FullMatrix<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
- FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
-
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
- M(i,j) += fe.shape_value (i, lhs_quadrature.point(q)) *
- fe.shape_value (j, lhs_quadrature.point(q)) *
- lhs_quadrature.weight(q);
-
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
- Q(i,q) += fe.shape_value (i, rhs_quadrature.point(q)) *
- rhs_quadrature.weight(q);
-
- // then invert M
- FullMatrix<double> 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 <int dim, int spacedim>
- void
- compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
- const Quadrature<dim> &quadrature,
- FullMatrix<double> &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<quadrature.size(); ++q)
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- I_q(q,i) = fe.shape_value (i, quadrature.point(q));
- }
-
-
-
- template <int dim>
- void
- compute_projection_from_quadrature_points(
- const FullMatrix<double> &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<double> component_at_node(n_support_points);
- // component at the quadrature point
- Vector<double> 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 <n_support_points; ++nn)
- {
- vector_of_tensors_at_nodes[nn][ii] = component_at_node(nn);
- }
- }
- }
-
-
-
- template <int dim>
- void
- compute_projection_from_quadrature_points(
- const FullMatrix<double> &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<double> component_at_node(n_support_points);
- // component at the quadrature point
- Vector<double> 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 <n_support_points; ++nn)
- {
- (vector_of_tensors_at_nodes[nn])[row][column] = component_at_node(nn);
- }
- }
- }
-
-
-
- template <int dim, int spacedim>
- void
- compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
- const Quadrature<dim-1> &lhs_quadrature,
- const Quadrature<dim-1> &rhs_quadrature,
- const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
- const unsigned int face,
- FullMatrix<double> &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<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
- FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
-
- {
- // need an FEFaceValues object to evaluate shape function
- // values on the specified face.
- FEFaceValues <dim> 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.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
- M(i,j) += fe_face_values.shape_value (i, q) *
- fe_face_values.shape_value (j, q) *
- lhs_quadrature.weight(q);
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- {
- M(i,i) = (M(i,i) == 0 ? 1 : M(i,i));
- }
- }
-
- {
- FEFaceValues <dim> 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<fe.dofs_per_cell; ++i)
- for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
- Q(i,q) += fe_face_values.shape_value (i, q) *
- rhs_quadrature.weight(q);
- }
- // then invert M
- FullMatrix<double> 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 <int dim>
- void
- hierarchic_to_lexicographic_numbering (unsigned int degree, std::vector<unsigned int> &h2l)
- {
- // number of support points in each direction
- const unsigned int n = degree+1;
-
- const unsigned int dofs_per_cell = Utilities::fixed_power<dim>(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<dofs_per_cell; ++i)
- h2l[i] = i-1;
-
- break;
- }
-
- case 2:
- {
- unsigned int next_index = 0;
- // first the four vertices
- h2l[next_index++] = 0;
- h2l[next_index++] = n-1;
- h2l[next_index++] = n*(n-1);
- h2l[next_index++] = n*n-1;
-
- // left line
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = (1+i)*n;
-
- // right line
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = (2+i)*n-1;
-
- // bottom line
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = 1+i;
-
- // top line
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = n*(n-1)+i+1;
-
- // inside quad
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- h2l[next_index++] = n*(i+1)+j+1;
-
- Assert (next_index == dofs_per_cell, ExcInternalError());
-
- break;
- }
-
- case 3:
- {
- unsigned int next_index = 0;
- // first the eight vertices
- h2l[next_index++] = 0; // 0
- h2l[next_index++] = ( 1)*degree; // 1
- h2l[next_index++] = ( n )*degree; // 2
- h2l[next_index++] = ( n+1)*degree; // 3
- h2l[next_index++] = (n*n )*degree; // 4
- h2l[next_index++] = (n*n +1)*degree; // 5
- h2l[next_index++] = (n*n+n )*degree; // 6
- h2l[next_index++] = (n*n+n+1)*degree; // 7
-
- // line 0
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = (i+1)*n;
- // line 1
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = n-1+(i+1)*n;
- // line 2
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = 1+i;
- // line 3
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = 1+i+n*(n-1);
-
- // line 4
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = (n-1)*n*n+(i+1)*n;
- // line 5
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = (n-1)*(n*n+1)+(i+1)*n;
- // line 6
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = n*n*(n-1)+i+1;
- // line 7
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = n*n*(n-1)+i+1+n*(n-1);
-
- // line 8
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = (i+1)*n*n;
- // line 9
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = n-1+(i+1)*n*n;
- // line 10
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = (i+1)*n*n+n*(n-1);
- // line 11
- for (unsigned int i=0; i<dofs_per_line; ++i)
- h2l[next_index++] = n-1+(i+1)*n*n+n*(n-1);
-
-
- // inside quads
- // face 0
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- h2l[next_index++] = (i+1)*n*n+n*(j+1);
- // face 1
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- h2l[next_index++] = (i+1)*n*n+n-1+n*(j+1);
- // face 2, note the orientation!
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- h2l[next_index++] = (j+1)*n*n+i+1;
- // face 3, note the orientation!
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- h2l[next_index++] = (j+1)*n*n+n*(n-1)+i+1;
- // face 4
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- h2l[next_index++] = n*(i+1)+j+1;
- // face 5
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- h2l[next_index++] = (n-1)*n*n+n*(i+1)+j+1;
-
- // inside hex
- for (unsigned int i=0; i<dofs_per_line; ++i)
- for (unsigned int j=0; j<dofs_per_line; ++j)
- for (unsigned int k=0; k<dofs_per_line; ++k)
- h2l[next_index++] = n*n*(i+1)+n*(j+1)+k+1;
-
- Assert (next_index == dofs_per_cell, ExcInternalError());
-
- break;
- }
-
- default:
- Assert (false, ExcNotImplemented());
- }
- }
-
-
-
- template <int dim>
- void
- hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &fe,
- std::vector<unsigned int> &h2l)
- {
- Assert (h2l.size() == fe.dofs_per_cell,
- ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell));
- hierarchic_to_lexicographic_numbering<dim> (fe.dofs_per_line+1, h2l);
- }
-
-
-
- template <int dim>
- std::vector<unsigned int>
- hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &fe)
- {
- Assert (fe.n_components() == 1, ExcInvalidFE());
- std::vector<unsigned int> h2l(fe.dofs_per_cell);
- hierarchic_to_lexicographic_numbering<dim> (fe.dofs_per_line+1, h2l);
- return (h2l);
- }
-
-
-
- template <int dim>
- void
- lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe,
- std::vector<unsigned int> &l2h)
- {
- l2h = lexicographic_to_hierarchic_numbering (fe);
- }
-
-
-
- template <int dim>
- std::vector<unsigned int>
- lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &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