void initialize (const std::vector<const FiniteElement<dim,spacedim>*> &fes,
const std::vector<unsigned int> &multiplicities);
- /**
- * Used by @p initialize.
- */
- void build_face_tables();
-
/**
* Used by @p initialize.
*/
const FiniteElement<dim,spacedim> &finite_element,
const bool do_tensor_product = true);
+ /**
+ * build...
+ * @param face_system_to_base_table
+ * @param face_system_to_component_table
+ * @param finite_element
+ * @param do_tensor_product
+ */
+ 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> &finite_element,
+ const bool do_tensor_product = true);
+
/**
* Parse the name of a finite element and generate a finite element object
* accordingly. The parser ignores space characters between words (things
}
-template <int dim, int spacedim>
-void
-FESystem<dim,spacedim>::build_face_tables()
-{
- // 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<this->n_base_elements(); ++base)
- for (unsigned int m=0; m<this->element_multiplicity(base);
- ++m, comp_start += base_element(base).n_components())
- for (unsigned int local_index = 0;
- local_index < 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
- = (base_element(base).dofs_per_vertex*vertex_number +
- local_index);
-
- const unsigned int face_index_in_base
- = (base_element(base).dofs_per_vertex*vertex_number +
- local_index);
-
- this->face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
- if (base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = 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
- = base_element(base).face_system_to_component_index(face_index_in_base).second;
- this->face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
- }
- else
- this->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 < this->n_base_elements(); ++base)
- for (unsigned int m=0; m<this->element_multiplicity(base);
- ++m, comp_start += base_element(base).n_components())
- for (unsigned int local_index = 0;
- local_index < base_element(base).dofs_per_line;
- ++local_index, ++total_index)
- {
- // do everything alike for this type of object
- const unsigned int index_in_base
- = (base_element(base).dofs_per_line*line_number +
- local_index +
- base_element(base).first_line_index);
-
- const unsigned int face_index_in_base
- = (base_element(base).first_face_line_index +
- base_element(base).dofs_per_line * line_number +
- local_index);
-
- this->face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
- if (base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = 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
- = base_element(base).face_system_to_component_index(face_index_in_base).second;
- this->face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
- }
- else
- this->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<this->n_base_elements(); ++base)
- for (unsigned int m=0; m<this->element_multiplicity(base);
- ++m, comp_start += base_element(base).n_components())
- for (unsigned int local_index = 0;
- local_index < base_element(base).dofs_per_quad;
- ++local_index, ++total_index)
- {
- // do everything alike for this type of object
- const unsigned int index_in_base
- = (base_element(base).dofs_per_quad*quad_number +
- local_index +
- base_element(base).first_quad_index);
-
- const unsigned int face_index_in_base
- = (base_element(base).first_face_quad_index +
- base_element(base).dofs_per_quad * quad_number +
- local_index);
-
- this->face_system_to_base_table[total_index]
- = std::make_pair (std::make_pair(base,m), face_index_in_base);
-
- if (base_element(base).is_primitive(index_in_base))
- {
- const unsigned int comp_in_base
- = 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
- = base_element(base).face_system_to_component_index(face_index_in_base).second;
- this->face_system_to_component_table[total_index]
- = std::make_pair (comp, face_index_in_comp);
- }
- else
- this->face_system_to_component_table[total_index] = non_primitive_index;
- }
- }
- Assert (total_index == this->dofs_per_face, ExcInternalError());
- Assert (total_index == this->face_system_to_component_table.size(),
- ExcInternalError());
- Assert (total_index == this->face_system_to_base_table.size(),
- ExcInternalError());
-}
-
-
-
template <int dim, int spacedim>
void FESystem<dim,spacedim>::build_interface_constraints ()
{
this->component_to_base_table,
*this);
+ FETools::build_face_tables(this->face_system_to_base_table,
+ this->face_system_to_component_table,
+ *this);
+
}
- build_face_tables();
// restriction and prolongation matrices are build on demand
}
}
+ 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())
+ 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())
+ 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())
+ 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> *
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > &component_to_base_table,
const FiniteElement<deal_II_dimension,deal_II_space_dimension> &fe,
const bool do_tensor_product);
+
+ template
+ 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<deal_II_dimension,deal_II_space_dimension> &fe,
+ const bool do_tensor_product);
template
void compute_block_renumbering (