};
#endif
- std::vector<unsigned int> dofs_on_this_face;
+ std::vector<types::global_dof_index> dofs_on_this_face;
dofs_on_this_face.reserve (max_dofs_per_face(dof));
- // loop over all faces to check
- // whether they are at a
- // boundary. note that we need not
- // take special care of single
- // lines (using
- // @p{cell->has_boundary_lines}),
- // since we do not support
- // boundaries of dimension dim-2,
- // and so every boundary line is
- // also part of a boundary face.
+ // loop over all faces to check whether they are at a boundary. note
+ // that we need not take special care of single lines (using
+ // @p{cell->has_boundary_lines}), since we do not support boundaries of
+ // dimension dim-2, and so every boundary line is also part of a
+ // boundary face.
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
for (; cell!=endc; ++cell)
cell = cell->child(direction);
const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex;
- std::vector<unsigned int> boundary_dof_boundary_indices (dofs_per_vertex);
+ std::vector<types::global_dof_index> boundary_dof_boundary_indices (dofs_per_vertex);
- // next get boundary mapped dof
- // indices of boundary dofs
+ // next get boundary mapped dof indices of boundary dofs
for (unsigned int i=0; i<dofs_per_vertex; ++i)
boundary_dof_boundary_indices[i]
= dof_to_boundary_mapping[cell->vertex_dof_index(direction,i)];
{
namespace
{
- inline
- bool
+ inline bool
check_master_dof_list (const FullMatrix<double> &face_interpolation_matrix,
- const std::vector<unsigned int> &master_dof_list)
+ const std::vector<types::global_dof_index> &master_dof_list)
{
const unsigned int N = master_dof_list.size();
(fe2.dofs_per_quad <= fe1.dofs_per_quad),
ExcInternalError());
- // the idea here is to designate as
- // many DoFs in fe1 per object
- // (vertex, line, quad) as master as
- // there are such dofs in fe2
- // (indices are int, because we want
- // to avoid the 'unsigned int < 0 is
- // always false warning for the cases
- // at the bottom in 1d and 2d)
+ // the idea here is to designate as many DoFs in fe1 per object
+ // (vertex, line, quad) as master as there are such dofs in fe2
+ // (indices are int, because we want to avoid the 'unsigned int < 0
+ // is always false warning for the cases at the bottom in 1d and
+ // 2d)
//
- // as mentioned in the paper, it is
- // not always easy to find a set of
- // master dofs that produces an
- // invertible matrix. to this end, we
- // check in each step whether the
- // matrix is still invertible and
- // simply discard this dof if the
- // matrix is not invertible anymore.
+ // as mentioned in the paper, it is not always easy to find a set
+ // of master dofs that produces an invertible matrix. to this end,
+ // we check in each step whether the matrix is still invertible and
+ // simply discard this dof if the matrix is not invertible anymore.
//
- // the cases where we did have
- // trouble in the past were with
- // adding more quad dofs when Q3 and
- // Q4 elements meet at a refined face
- // in 3d (see the hp/crash_12 test
- // that tests that we can do exactly
- // this, and failed before we had
- // code to compensate for this
- // case). the other case are system
- // elements: if we have say a Q1Q2 vs
- // a Q2Q3 element, then we can't just
- // take all master dofs on a line
- // from a single base element, since
- // the shape functions of that base
- // element are independent of that of
- // the other one. this latter case
- // shows up when running
- // hp/hp_constraints_q_system_06
+ // the cases where we did have trouble in the past were with adding
+ // more quad dofs when Q3 and Q4 elements meet at a refined face in
+ // 3d (see the hp/crash_12 test that tests that we can do exactly
+ // this, and failed before we had code to compensate for this
+ // case). the other case are system elements: if we have say a Q1Q2
+ // vs a Q2Q3 element, then we can't just take all master dofs on a
+ // line from a single base element, since the shape functions of
+ // that base element are independent of that of the other one. this
+ // latter case shows up when running hp/hp_constraints_q_system_06
+
- std::vector<unsigned int> master_dof_list;
+ std::vector<types::global_dof_index> master_dof_list;
unsigned int index = 0;
for (int v=0;
v<static_cast<signed int>(GeometryInfo<dim>::vertices_per_face);
AssertDimension (index, fe1.dofs_per_face);
AssertDimension (master_dof_list.size(), fe2.dofs_per_face);
- // finally copy the list into the
- // mask
+ // finally copy the list into the mask
std::fill (master_dof_mask.begin(), master_dof_mask.end(), false);
- for (std::vector<unsigned int>::const_iterator i=master_dof_list.begin();
+ for (std::vector<types::global_dof_index>::const_iterator i=master_dof_list.begin();
i!=master_dof_list.end(); ++i)
master_dof_mask[*i] = true;
}
/**
- * Copy constraints into a constraint
- * matrix object.
+ * Copy constraints into a constraint matrix object.
*
- * This function removes zero
- * constraints and those, which
- * constrain a DoF which was
- * already eliminated in one of
- * the previous steps of the hp
- * hanging node procedure.
+ * This function removes zero constraints and those, which constrain
+ * a DoF which was already eliminated in one of the previous steps of
+ * the hp hanging node procedure.
*
- * It also suppresses very small
- * entries in the constraint matrix to
- * avoid making the sparsity pattern
- * fuller than necessary.
+ * It also suppresses very small entries in the constraint matrix to
+ * avoid making the sparsity pattern fuller than necessary.
*/
void
- filter_constraints (const std::vector<unsigned int> &master_dofs,
- const std::vector<unsigned int> &slave_dofs,
+ filter_constraints (const std::vector<types::global_dof_index> &master_dofs,
+ const std::vector<types::global_dof_index> &slave_dofs,
const FullMatrix<double> &face_constraints,
ConstraintMatrix &constraints)
{
void
- make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> & /*dof_handler*/,
- ConstraintMatrix & /*constraints*/)
+ make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/,
+ ConstraintMatrix &/*constraints*/)
{
- // we may have to compute
- // constraints for
- // vertices. gotta think about
+ // we may have to compute constraints for vertices. gotta think about
// that a bit more
- //TODO[WB]: think about what to do here...
+
+ //TODO[WB]: think about what to do here...
}
void
- make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> & /*dof_handler*/,
- ConstraintMatrix & /*constraints*/,
+ make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/,
+ ConstraintMatrix &/*constraints*/,
dealii::internal::int2type<1>)
{
- // we may have to compute
- // constraints for
- // vertices. gotta think about
+ // we may have to compute constraints for vertices. gotta think about
// that a bit more
- //TODO[WB]: think about what to do here...
+
+ //TODO[WB]: think about what to do here...
}
const unsigned int spacedim = DH::space_dimension;
- std::vector<unsigned int> dofs_on_mother;
- std::vector<unsigned int> dofs_on_children;
+ std::vector<types::global_dof_index> dofs_on_mother;
+ std::vector<types::global_dof_index> dofs_on_children;
- // loop over all lines; only on
- // lines there can be constraints.
- // We do so by looping over all
- // active cells and checking
- // whether any of the faces are
- // refined which can only be from
- // the neighboring cell because
- // this one is active. In that
- // case, the face is subject to
+ // loop over all lines; only on lines there can be constraints. We do
+ // so by looping over all active cells and checking whether any of
+ // the faces are refined which can only be from the neighboring cell
+ // because this one is active. In that case, the face is subject to
// constraints
//
- // note that even though we may
- // visit a face twice if the
- // neighboring cells are equally
- // refined, we can only visit each
- // face with hanging nodes once
+ // note that even though we may visit a face twice if the neighboring
+ // cells are equally refined, we can only visit each face with
+ // hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
const unsigned int dim = 3;
- std::vector<unsigned int> dofs_on_mother;
- std::vector<unsigned int> dofs_on_children;
+ std::vector<types::global_dof_index> dofs_on_mother;
+ std::vector<types::global_dof_index> dofs_on_children;
- // loop over all quads; only on
- // quads there can be constraints.
- // We do so by looping over all
- // active cells and checking
- // whether any of the faces are
- // refined which can only be from
- // the neighboring cell because
- // this one is active. In that
- // case, the face is subject to
+ // loop over all quads; only on quads there can be constraints. We do
+ // so by looping over all active cells and checking whether any of
+ // the faces are refined which can only be from the neighboring cell
+ // because this one is active. In that case, the face is subject to
// constraints
//
- // note that even though we may
- // visit a face twice if the
- // neighboring cells are equally
- // refined, we can only visit each
- // face with hanging nodes once
+ // note that even though we may visit a face twice if the neighboring
+ // cells are equally refined, we can only visit each face with
+ // hanging nodes once
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
const unsigned int spacedim = DH::space_dimension;
- // a matrix to be used for
- // constraints below. declared
- // here and simply resized down
- // below to avoid permanent
- // re-allocation of memory
+ // a matrix to be used for constraints below. declared here and
+ // simply resized down below to avoid permanent re-allocation of
+ // memory
FullMatrix<double> constraint_matrix;
- // similarly have arrays that
- // will hold master and slave
- // dof numbers, as well as a
- // scratch array needed for the
- // complicated case below
+ // similarly have arrays that will hold master and slave dof numbers,
+ // as well as a scratch array needed for the complicated case below
- std::vector<unsigned int> master_dofs;
- std::vector<unsigned int> slave_dofs;
- std::vector<unsigned int> scratch_dofs;
+ std::vector<types::global_dof_index> master_dofs;
+ std::vector<types::global_dof_index> slave_dofs;
+ std::vector<types::global_dof_index> scratch_dofs;
- // caches for the face and
- // subface interpolation
- // matrices between different
- // (or the same) finite
- // elements. we compute them
- // only once, namely the first
- // time they are needed, and
- // then just reuse them
+ // caches for the face and subface interpolation matrices between
+ // different (or the same) finite elements. we compute them only
+ // once, namely the first time they are needed, and then just reuse
+ // them
Table<2,std_cxx1x::shared_ptr<FullMatrix<double> > >
face_interpolation_matrices (n_finite_elements (dof_handler),
n_finite_elements (dof_handler));
ExcInternalError());
}
- // then loop over all cells and do
- // the work
+ // then loop over all cells and do the work
- std::vector<unsigned int> indices;
+ std::vector<types::global_dof_index> indices;
for (typename DH::active_cell_iterator c=dof.begin_active();
c!=dof.end(); ++ c)
if (c->is_locally_owned())
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
local_selected_dofs[i] = component_mask[local_component_asssociation[i]];
- // then loop over all cells and do
- // work
+ // then loop over all cells and do work
- std::vector<unsigned int> indices(fe.dofs_per_cell);
- typename MGDoFHandler<dim,spacedim>::cell_iterator c;
+ std::vector<types::global_dof_index> indices(fe.dofs_per_cell);
+ typename DH::cell_iterator c;
for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
{
c->get_mg_dof_indices(indices);
(component_mask.n_selected_components(n_components(dof_handler)) !=
n_components(dof_handler)));
- std::vector<unsigned int> face_dof_indices;
+ std::vector<types::global_dof_index> face_dof_indices;
face_dof_indices.reserve (max_dofs_per_face(dof_handler));
- // now loop over all cells and
- // check whether their faces are at
- // the boundary. note that we need
- // not take special care of single
- // lines being at the boundary
- // (using
- // @p{cell->has_boundary_lines}),
- // since we do not support
- // boundaries of dimension dim-2,
- // and so every isolated boundary
- // line is also part of a boundary
- // face which we will be visiting
- // sooner or later
+ // now loop over all cells and check whether their faces are at the
+ // boundary. note that we need not take special care of single lines
+ // being at the boundary (using @p{cell->has_boundary_lines}), since we
+ // do not support boundaries of dimension dim-2, and so every isolated
+ // boundary line is also part of a boundary face which we will be
+ // visiting sooner or later
for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end(); ++cell)
- // only work on cells that are either locally owned or at
- // least ghost cells
+
+ // only work on cells that are either locally owned or at least ghost
+ // cells
if (cell->is_artificial() == false)
for (unsigned int face=0;
face<GeometryInfo<DH::dimension>::faces_per_cell; ++face)
const bool check_vector_component
= (component_mask.represents_the_all_selected_mask() == false);
- // clear and reset array by default
- // values
+ // clear and reset array by default values
selected_dofs.clear ();
selected_dofs.resize (dof_handler.n_dofs(), false);
- std::vector<unsigned int> cell_dof_indices;
+ std::vector<types::global_dof_index> cell_dof_indices;
cell_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // now loop over all cells and
- // check whether their faces are at
- // the boundary. note that we need
- // not take special care of single
- // lines being at the boundary
- // (using
- // @p{cell->has_boundary_lines}),
- // since we do not support
- // boundaries of dimension dim-2,
- // and so every isolated boundary
- // line is also part of a boundary
- // face which we will be visiting
- // sooner or later
+ // now loop over all cells and check whether their faces are at the
+ // boundary. note that we need not take special care of single lines
+ // being at the boundary (using @p{cell->has_boundary_lines}), since we
+ // do not support boundaries of dimension dim-2, and so every isolated
+ // boundary line is also part of a boundary face which we will be
+ // visiting sooner or later
for (typename DH::active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end(); ++cell)
for (unsigned int face=0;
// preset all values by false
std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
- std::vector<unsigned int> local_dof_indices;
+ std::vector<types::global_dof_index> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // this function is similar to the
- // make_sparsity_pattern function,
- // see there for more information
+ // this function is similar to the make_sparsity_pattern function, see
+ // there for more information
typename DH::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
// collect all the locally owned dofs
dof_set = dof_handler.locally_owned_dofs();
- // add the DoF on the adjacent ghost cells
- // to the IndexSet, cache them in a
- // set. need to check each dof manually
- // because we can't be sure that the dof
- // range of locally_owned_dofs is really
- // contiguous.
+ // add the DoF on the adjacent ghost cells to the IndexSet, cache them
+ // in a set. need to check each dof manually because we can't be sure
+ // that the dof range of locally_owned_dofs is really contiguous.
- std::vector<unsigned int> dof_indices;
- std::set<unsigned int> global_dof_indices;
+ std::vector<types::global_dof_index> dof_indices;
+ std::set<types::global_dof_index> global_dof_indices;
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
// collect all the locally owned dofs
dof_set = dof_handler.locally_owned_dofs();
- // add the DoF on the adjacent ghost cells
- // to the IndexSet, cache them in a
- // set. need to check each dof manually
- // because we can't be sure that the dof
- // range of locally_owned_dofs is really
- // contiguous.
+ // add the DoF on the adjacent ghost cells to the IndexSet, cache them
+ // in a set. need to check each dof manually because we can't be sure
+ // that the dof range of locally_owned_dofs is really contiguous.
- std::vector<unsigned int> dof_indices;
- std::set<unsigned int> global_dof_indices;
+ std::vector<types::global_dof_index> dof_indices;
+ std::set<types::global_dof_index> global_dof_indices;
typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
std::fill_n (subdomain_association.begin(), dof_handler.n_dofs(),
types::invalid_subdomain_id);
- std::vector<unsigned int> local_dof_indices;
+ std::vector<types::global_dof_index> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // pseudo-randomly assign variables
- // which lie on the interface
- // between subdomains to each of
- // the two or more
+ // pseudo-randomly assign variables which lie on the interface between
+ // subdomains to each of the two or more
bool coin_flip = true;
- // loop over all cells and record
- // which subdomain a DoF belongs
- // to. toss a coin in case it is on
- // an interface
+ // loop over all cells and record which subdomain a DoF belongs to.
+ // toss a coin in case it is on an interface
typename DH::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
IndexSet index_set (dof_handler.n_dofs());
- std::vector<unsigned int> local_dof_indices;
+ std::vector<types::global_dof_index> local_dof_indices;
local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
- // first generate an unsorted list of all
- // indices which we fill from the back. could
- // also insert them directly into the
- // IndexSet, but that inserts indices in the
- // middle, which is an O(n^2) algorithm and
- // hence too expensive. Could also use
- // std::set, but that is in general more
- // expensive than a vector
+ // first generate an unsorted list of all indices which we fill from
+ // the back. could also insert them directly into the IndexSet, but
+ // that inserts indices in the middle, which is an O(n^2) algorithm and
+ // hence too expensive. Could also use std::set, but that is in general
+ // more expensive than a vector
- std::vector<unsigned int> subdomain_indices;
+ std::vector<types::global_dof_index> subdomain_indices;
typename DH::active_cell_iterator
cell = dof_handler.begin_active(),
void
resolve_components (const hp::FECollection<dim,spacedim> &fe_collection,
const std::vector<unsigned char> &dofs_by_component,
- const std::vector<unsigned int> &target_component,
+ const std::vector<unsigned int> &target_component,
const bool only_once,
- std::vector<unsigned int> &dofs_per_component,
+ std::vector<types::global_dof_index> &dofs_per_component,
unsigned int &component)
{
- // assert that all elements in the collection have the same
- // structure (base elements and multiplicity, components per base
- // element) and then simply call the function above
+ // assert that all elements in the collection have the same structure
+ // (base elements and multiplicity, components per base element) and
+ // then simply call the function above
for (unsigned int fe=1; fe<fe_collection.size(); ++fe)
{
Assert (fe_collection[fe].n_components() == fe_collection[0].n_components(),
unsigned int n_parameters_on_fine_grid=0;
if (true)
{
- // have a flag for each dof on
- // the fine grid and set it
- // to true if this is an
- // interesting dof. finally count
- // how many true's there
+ // have a flag for each dof on the fine grid and set it to true
+ // if this is an interesting dof. finally count how many true's
+ // there
std::vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
- std::vector<unsigned int> local_dof_indices (fine_fe.dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices (fine_fe.dofs_per_cell);
for (typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
cell=fine_grid.begin_active();
template <class InputVector, typename number>
void FEValuesBase<dim,spacedim>::get_function_values (
const InputVector &fe_function,
- const VectorSlice<const std::vector<unsigned int> > &indices,
+ const VectorSlice<const std::vector<types::global_dof_index> > &indices,
std::vector<Vector<number> > &values) const
{
- // One value per quadrature point
- Assert (n_quadrature_points == values.size(),
- ExcDimensionMismatch(values.size(), n_quadrature_points));
-
- const unsigned int n_components = fe->n_components();
-
- // Size of indices must be a
- // multiple of dofs_per_cell such
- // that an integer number of
- // function values is generated in
- // each point.
+ // Size of indices must be a multiple of dofs_per_cell such that an integer
+ // number of function values is generated in each point.
Assert (indices.size() % dofs_per_cell == 0,
ExcNotMultiple(indices.size(), dofs_per_cell));
-
- // The number of components of the
- // result may be a multiple of the
- // number of components of the
- // finite element
- const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
-
- for (unsigned int i=0; i<values.size(); ++i)
- Assert (values[i].size() == result_components,
- ExcDimensionMismatch(values[i].size(), result_components));
-
- // If the result has more
- // components than the finite
- // element, we need this number for
- // loops filling all components
- const unsigned int component_multiple = result_components / n_components;
-
Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
- // initialize with zero
- for (unsigned int i=0; i<values.size(); ++i)
- std::fill_n (values[i].begin(), values[i].size(), 0);
-
- // add up contributions of trial
- // functions. now check whether the
- // shape function is primitive or
- // not. if it is, then set its only
- // non-zero component, otherwise
- // loop over components
- for (unsigned int mc = 0; mc < component_multiple; ++mc)
- for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
- {
- const double value = get_vector_element (fe_function, indices[shape_func+mc*dofs_per_cell]);
- if (value == 0.)
- continue;
-
- if (fe->is_primitive(shape_func))
- {
- const unsigned int comp = fe->system_to_component_index(shape_func).first
- + mc * n_components;
- const unsigned int
- row = this->shape_function_to_row_table[shape_func * fe->n_components() + comp];
-
- const double *shape_value_ptr = &this->shape_values(row, 0);
- for (unsigned int point=0; point<n_quadrature_points; ++point)
- values[point](comp) += value * *shape_value_ptr++;
- }
- else
- for (unsigned int c=0; c<n_components; ++c)
- {
- if (fe->get_nonzero_components(shape_func)[c] == false)
- continue;
-
- const unsigned int
- row = this->shape_function_to_row_table[shape_func * fe->n_components() + c];
-
- const double *shape_value_ptr = &this->shape_values(row, 0);
- const unsigned int comp = c + mc * n_components;
-
- for (unsigned int point=0; point<n_quadrature_points; ++point)
- values[point](comp) += value * *shape_value_ptr++;
- }
- }
+ VectorSlice<std::vector<Vector<number> > > val(values);
+ if (indices.size() <= 100)
+ {
+ double dof_values[100];
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ dof_values[i] = get_vector_element (fe_function, indices[i]);
+ internal::do_function_values(&dof_values[0], this->shape_values, *fe,
+ this->shape_function_to_row_table, val,
+ false, indices.size()/dofs_per_cell);
+ }
+ else
+ {
+ Vector<double> dof_values(100);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ dof_values[i] = get_vector_element (fe_function, indices[i]);
+ internal::do_function_values(dof_values.begin(), this->shape_values, *fe,
+ this->shape_function_to_row_table, val,
+ false, indices.size()/dofs_per_cell);
+ }
}
template <class InputVector>
void FEValuesBase<dim,spacedim>::get_function_hessians (
const InputVector &fe_function,
- const VectorSlice<const std::vector<unsigned int> > &indices,
+ const VectorSlice<const std::vector<types::global_dof_index> > &indices,
std::vector<Tensor<2,spacedim> > &hessians) const
{
- Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField());
- // This function fills a single
- // component only
- Assert (fe->n_components() == 1,
- ExcDimensionMismatch(fe->n_components(), 1));
- // One index for each dof
- Assert (indices.size() == dofs_per_cell,
- ExcDimensionMismatch(indices.size(), dofs_per_cell));
- // This vector has one entry for
- // each quadrature point
- Assert (hessians.size() == n_quadrature_points,
- ExcDimensionMismatch(hessians.size(), n_quadrature_points));
-
- // initialize with zero
- std::fill_n (hessians.begin(), n_quadrature_points, Tensor<2,spacedim>());
-
- // add up contributions of trial
- // functions. note that here we
- // deal with scalar finite
- // elements, so no need to check
- // for non-primitivity of shape
- // functions
- for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+ Assert (this->update_flags & update_second_derivatives,
+ ExcAccessToUninitializedField());
+ AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler());
+ AssertDimension (indices.size(), dofs_per_cell);
+ if (dofs_per_cell <= 100)
+ {
+ double dof_values[100];
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ dof_values[i] = get_vector_element (fe_function, indices[i]);
+ internal::do_function_derivatives(&dof_values[0], this->shape_hessians,
+ hessians);
+ }
+ else
{
- const double value = get_vector_element (fe_function, indices[shape_func]);
- if (value == 0.)
- continue;
-
- const Tensor<2,spacedim> *shape_hessians_ptr
- = &this->shape_hessians[shape_func][0];
- for (unsigned int point=0; point<n_quadrature_points; ++point)
- hessians[point] += value * *shape_hessians_ptr++;
+ Vector<double> dof_values(dofs_per_cell);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ dof_values[i] = get_vector_element (fe_function, indices[i]);
+ internal::do_function_derivatives(dof_values.begin(), this->shape_hessians,
+ hessians);
}
}
template <class InputVector, typename number>
void FEValuesBase<dim,spacedim>::get_function_laplacians (
const InputVector &fe_function,
- const VectorSlice<const std::vector<unsigned int> > &indices,
+ const VectorSlice<const std::vector<types::global_dof_index> > &indices,
std::vector<Vector<number> > &laplacians) const
{
- // One value per quadrature point
- Assert (n_quadrature_points == laplacians.size(),
- ExcDimensionMismatch(laplacians.size(), n_quadrature_points));
-
- const unsigned int n_components = fe->n_components();
-
- // Size of indices must be a
- // multiple of dofs_per_cell such
- // that an integer number of
- // function values is generated in
- // each point.
+ // Size of indices must be a multiple of dofs_per_cell such that an integer
+ // number of function values is generated in each point.
Assert (indices.size() % dofs_per_cell == 0,
ExcNotMultiple(indices.size(), dofs_per_cell));
-
- // The number of components of the
- // result may be a multiple of the
- // number of components of the
- // finite element
- const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
-
- for (unsigned int i=0; i<laplacians.size(); ++i)
- Assert (laplacians[i].size() == result_components,
- ExcDimensionMismatch(laplacians[i].size(), result_components));
-
- // If the result has more
- // components than the finite
- // element, we need this number for
- // loops filling all components
- const unsigned int component_multiple = result_components / n_components;
-
Assert (this->update_flags & update_hessians, ExcAccessToUninitializedField());
-
- // initialize with zero
- for (unsigned int i=0; i<laplacians.size(); ++i)
- std::fill_n (laplacians[i].begin(), laplacians[i].size(), 0);
-
- // add up contributions of trial
- // functions. now check whether the
- // shape function is primitive or
- // not. if it is, then set its only
- // non-zero component, otherwise
- // loop over components
- for (unsigned int mc = 0; mc < component_multiple; ++mc)
- for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
- {
- const double value = get_vector_element (fe_function, indices[shape_func+mc*dofs_per_cell]);
- if (value == 0.)
- continue;
-
- if (fe->is_primitive(shape_func))
- {
- const unsigned int comp = fe->system_to_component_index(shape_func).first
- + mc * n_components;
- const unsigned int
- row = this->shape_function_to_row_table[shape_func * fe->n_components() + comp];
-
- const Tensor<2,spacedim> *shape_hessian_ptr
- = &this->shape_hessians[row][0];
- for (unsigned int point=0; point<n_quadrature_points; ++point)
- laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
- }
- else
- for (unsigned int c=0; c<n_components; ++c)
- {
- if (fe->get_nonzero_components(shape_func)[c] == false)
- continue;
-
- const unsigned int
- row = this->shape_function_to_row_table[shape_func * fe->n_components() + c];
-
- const Tensor<2,spacedim> *shape_hessian_ptr
- = &this->shape_hessians[row][0];
- const unsigned int comp = c + mc * n_components;
-
- for (unsigned int point=0; point<n_quadrature_points; ++point)
- laplacians[point](comp) += value * trace(*shape_hessian_ptr++);
- }
- }
+ if (indices.size() <= 100)
+ {
+ double dof_values[100];
+ for (unsigned int i=0; i<indices.size(); ++i)
+ dof_values[i] = get_vector_element (fe_function, indices[i]);
+ internal::do_function_laplacians(&dof_values[0], this->shape_hessians,
+ *fe, this->shape_function_to_row_table,
+ laplacians, false,
+ indices.size()/dofs_per_cell);
+ }
+ else
+ {
+ Vector<double> dof_values(indices.size());
+ for (unsigned int i=0; i<indices.size(); ++i)
+ dof_values[i] = get_vector_element (fe_function, indices[i]);
+ internal::do_function_laplacians(dof_values.begin(),this->shape_hessians,
+ *fe, this->shape_function_to_row_table,
+ laplacians, false,
+ indices.size()/dofs_per_cell);
+ }
}
projects / dealii-svn.git / commitdiff