*/
bool operator == (const FiniteElementBase<dim> &) const;
- /**
- * Given a vector component and
- * an index of a shape function
- * within the shape functions
- * corresponding to this vector
- * component, return the index of
- * this shape function within the
- * shape functions of this
- * element. If this is a scalar
- * element, then the given
- * component may only be zero,
- * and the given component index
- * is also the return value.
- *
- * If the finite element is
- * vector-valued and has
- * non-primitive shape functions,
- * i.e. some of its shape
- * functions are non-zero in more
- * than just one vector
- * component, then this function
- * cannot be used since shape
- * functions are no more
- * associated with individual
- * vector components, and an
- * exception of type
- * @p{ExcFENotPrimitive} is
- * thrown.
- */
- unsigned int component_to_system_index (const unsigned int component,
- const unsigned int component_index) const;
-
- /**
- * Same as above, but compute the
- * data from the index of a shape
- * function on a face.
- */
- unsigned int face_component_to_system_index (const unsigned int component,
- const unsigned int component_index) const;
/**
* Compute vector component and
std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
face_system_to_base_table;
- /**
- * Map between component and
- * linear dofs: For each pair of
- * vector component and index
- * within this component, store
- * the global dof number in the
- * composed element. If the
- * element is scalar, then the
- * outer (component) index can
- * only be zero, and the inner
- * index is equal to the stored
- * value.
- *
- * If the element is not
- * primitive, i.e. there are
- * shape functions that are
- * non-zero in more than one
- * vector-component, then this
- * function is obviously useless,
- * and all entries will be
- * invalid.
- */
- std::vector< std::vector<unsigned int> > component_to_system_table;
-
- /**
- * Map between component and
- * linear dofs on a face. Same
- * applies as above.
- */
- std::vector< std::vector<unsigned int> > face_component_to_system_table;
-
/**
* The base element establishing
* a component.
};
-template <int dim>
-inline
-unsigned int
-FiniteElementBase<dim>::component_to_system_index (const unsigned int component,
- const unsigned int component_index) const
-{
- Assert(component<this->n_components(),
- ExcIndexRange(component, 0, this->n_components()));
- Assert(component_index<component_to_system_table[component].size(),
- ExcIndexRange(component_index, 0,
- component_to_system_table[component].size()));
- Assert (is_primitive(),
- typename FiniteElementBase<dim>::ExcFENotPrimitive());
-
- return component_to_system_table[component][component_index];
-}
-
template <int dim>
inline
}
-template <int dim>
-inline
-unsigned int
-FiniteElementBase<dim>::face_component_to_system_index (const unsigned int component,
- const unsigned int component_index) const
-{
- Assert(component<this->n_components(),
- ExcIndexRange(component, 0, this->n_components()));
- Assert(component_index<face_component_to_system_table[component].size(),
- ExcIndexRange(component_index, 0,
- face_component_to_system_table[component].size()));
- Assert (is_primitive(),
- typename FiniteElementBase<dim>::ExcFENotPrimitive());
-
- return face_component_to_system_table[component][component_index];
-}
template <int dim>
Assert (face_quadrature != 0, ExcInternalError());
}
-
- for (unsigned int base_no=0, comp=0; base_no<n_base_elements(); ++base_no)
+ // let base elements update the
+ // necessary data
+ for (unsigned int base_no=0; base_no<n_base_elements(); ++base_no)
{
- const FiniteElement<dim> &base_fe=base_element(base_no);
+ const FiniteElement<dim> &
+ base_fe = base_element(base_no);
typename FiniteElementBase<dim>::InternalDataBase &
base_fe_data = fe_data.get_fe_data(base_no);
-
- // Make sure that in the
- // case of fill_fe_values
- // the data is only copied
- // from base_data to data
- // if base_data is
- // changed. therefore use
- // fe_fe_data.current_update_flags()
-
- // for the case of
- // fill_fe_(sub)face_values
- // the data needs to be
- // copied from base_data to
- // data on each face,
- // therefore use
- // base_fe_data.update_flags.
-
- // Store these flags into
- // base_flags before
- // calling
- // base_fe.fill_fe_([sub]face_)values
- // as the latter changes
- // the return value of
- // base_fe_data.current_update_flags()
- const UpdateFlags base_flags(dim_1==dim ?
- base_fe_data.current_update_flags() :
- base_fe_data.update_flags);
-
- FEValuesData<dim> &base_data=fe_data.get_fe_values_data(base_no);
-
+ FEValuesData<dim> &
+ base_data = fe_data.get_fe_values_data(base_no);
+
if (face_no==invalid_face_number)
base_fe.fill_fe_values(mapping, cell,
*cell_quadrature, mapping_data, base_fe_data, base_data);
else
base_fe.fill_fe_subface_values(mapping, cell, face_no, sub_no,
*face_quadrature, mapping_data, base_fe_data, base_data);
+ };
- for (unsigned int m=0; m<element_multiplicity(base_no); ++m, ++comp)
- for (unsigned int point=0; point<quadrature.n_quadrature_points; ++point)
- for (unsigned int k=0; k<base_fe.dofs_per_cell; ++k)
- {
- const unsigned int system_index=this->component_to_system_index(comp,k);
- if (base_flags & update_values)
- data.shape_values(system_index, point)=
- base_data.shape_values(k,point);
-
- if (base_flags & update_gradients)
- data.shape_gradients[system_index][point]=
- base_data.shape_gradients[k][point];
-
- if (base_flags & update_second_derivatives)
- data.shape_2nd_derivatives[system_index][point]=
- base_data.shape_2nd_derivatives[k][point];
- }
- }
+ // now data has been generated,
+ // so copy it. we used to work
+ // by looping over all base
+ // elements, then over
+ // multiplicity, then over the
+ // shape functions from that
+ // base element, but that
+ // requires that we can infer
+ // the global number of a shape
+ // function from its number in
+ // the base element. for that
+ // we had the
+ // component_to_system_table.
+ //
+ // however, this does of course
+ // no longer work since we have
+ // non-primitive elements. so
+ // we go the other way round:
+ // loop over all shape
+ // functions of the composed
+ // element, and find from where
+ // to copy the data. to be able
+ // to cache some data, make
+ // things a little bit more
+ // complicated: loop over all
+ // base elements, then over all
+ // shape functions of the
+ // composed element, and only
+ // treat those shape functions
+ // that belong to a given base
+ // element
+ for (unsigned int base_no=0; base_no<n_base_elements(); ++base_no)
+ for (unsigned int system_index=0; system_index<this->dofs_per_cell;
+ ++system_index)
+ if (this->system_to_base_table[system_index].first.first == base_no)
+ {
+ const FiniteElement<dim> &
+ base_fe = base_element(base_no);
+ typename FiniteElementBase<dim>::InternalDataBase &
+ base_fe_data = fe_data.get_fe_data(base_no);
+ FEValuesData<dim> &
+ base_data = fe_data.get_fe_values_data(base_no);
+
+ const unsigned int n_q_points = quadrature.n_quadrature_points;
+
+ // Make sure that in
+ // the case of
+ // fill_fe_values the
+ // data is only copied
+ // from base_data to
+ // data if base_data is
+ // changed. therefore
+ // use
+ // fe_fe_data.current_update_flags()
+ //
+ // for the case of
+ // fill_fe_(sub)face_values
+ // the data needs to be
+ // copied from
+ // base_data to data on
+ // each face, therefore
+ // use
+ // base_fe_data.update_flags.
+ //
+ // Store these flags into
+ // base_flags before
+ // calling
+ // base_fe.fill_fe_([sub]face_)values
+ // as the latter changes
+ // the return value of
+ // base_fe_data.current_update_flags()
+ const UpdateFlags base_flags(dim_1==dim ?
+ base_fe_data.current_update_flags() :
+ base_fe_data.update_flags);
+
+ const unsigned int k = system_to_base_table[system_index].second;
+ Assert (k<base_fe.dofs_per_cell, ExcInternalError());
+
+ // now copy. if the
+ // shape function is
+ // primitive, then
+ // there is only one
+ // value to be copied,
+ // but for
+ // non-primitive
+ // elements, there
+ // might be more values
+ // to be copied
+ //
+ // so, find out from
+ // which index to take
+ // this one value, and
+ // to which index to
+ // put
+ unsigned int out_index = 0;
+ for (unsigned int i=0; i<system_index; ++i)
+ out_index += this->n_nonzero_components(i);
+ unsigned int in_index = 0;
+ for (unsigned int i=0; i<k; ++i)
+ in_index += base_fe.n_nonzero_components(i);
+
+ // then loop over the
+ // number of components
+ // to be copied
+ Assert (this->n_nonzero_components(system_index) ==
+ base_fe.n_nonzero_components(k),
+ ExcInternalError());
+ for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)
+ {
+ if (base_flags & update_values)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ data.shape_values[out_index+s][q] =
+ base_data.shape_values(in_index+s,q);
+
+ if (base_flags & update_gradients)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ data.shape_gradients[out_index+s][q]=
+ base_data.shape_gradients[in_index+s][q];
+
+ if (base_flags & update_second_derivatives)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ data.shape_2nd_derivatives[out_index+s][q]=
+ base_data.shape_2nd_derivatives[in_index+s][q];
+ };
+ };
if (fe_data.first_cell)
this->system_to_component_table[total_index] = non_primitive_index;
}
}
-
- // Initialize mapping from
- // components to linear
- // index. Fortunately, this is the
- // inverse of what we just did.
- //
- // note that we can only do this if
- // all sub-elements are primiive,
- // otherwise fill the table with
- // invalid values
- if (this->is_primitive())
- {
- std::vector<unsigned int> dofs_per_component (this->n_components(), 0);
- for (unsigned int sys=0; sys<this->dofs_per_cell; ++sys)
- ++dofs_per_component[this->system_to_component_index(sys).first];
- for (unsigned int component=0; component<this->n_components(); ++component)
- this->component_to_system_table[component].resize(dofs_per_component[component]);
-
- // then go the reverse way to fill the array
- for (unsigned int sys=0; sys<this->dofs_per_cell; ++sys)
- {
- const unsigned int
- comp = this->system_to_component_index(sys).first,
- index_in_comp = this->system_to_component_index(sys).second;
-
- Assert (comp < this->component_to_system_table.size(),
- ExcInternalError());
- Assert (index_in_comp < this->component_to_system_table[comp].size(),
- ExcInternalError());
- this->component_to_system_table[comp][index_in_comp] = sys;
- };
- }
- else
- // element is not primitive, so
- // fill elements with nonsense
- {
- std::vector< std::vector<unsigned int> > tmp;
- this->component_to_system_table.swap (tmp);
- };
};
ExcInternalError());
Assert (total_index == this->face_system_to_base_table.size(),
ExcInternalError());
-
- // finally, initialize reverse
- // mapping. same applies
- // w.r.t. primitivity as above:
- if (this->is_primitive())
- {
- std::vector<unsigned int> dofs_per_component (this->n_components(), 0);
- for (unsigned int sys=0; sys<this->dofs_per_face; ++sys)
- ++dofs_per_component[this->face_system_to_component_index(sys).first];
- for (unsigned int component=0; component<this->n_components(); ++component)
- this->face_component_to_system_table[component].resize(dofs_per_component[component]);
-
- // then go the reverse way to fill
- // the array
- for (unsigned int sys=0; sys<this->dofs_per_face; ++sys)
- {
- const unsigned int
- comp = this->face_system_to_component_index(sys).first,
- index_in_comp = this->face_system_to_component_index(sys).second;
-
- Assert (comp < this->face_component_to_system_table.size(),
- ExcInternalError());
- Assert (index_in_comp < this->face_component_to_system_table[comp].size(),
- ExcInternalError());
- this->face_component_to_system_table[comp][index_in_comp] = sys;
- };
- }
- else
- // element is not primitive, so
- // fill elements with nonsense
- {
- std::vector< std::vector<unsigned int> > tmp;
- this->face_component_to_system_table.swap (tmp);
- };
};