-
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::
-InternalDataBase::initialize_2nd (const FiniteElement<dim,spacedim> *element,
- const Mapping<dim,spacedim> &mapping,
- const Quadrature<dim> &quadrature)
-{
- // if we shall compute second
- // derivatives, then we do so by
- // finite differencing the
- // gradients. that we do by
- // evaluating the gradients of
- // shape values at points shifted
- // star-like a little in each
- // coordinate direction around each
- // quadrature point.
- //
- // therefore generate 2*dim (the
- // number of evaluation points)
- // FEValues objects with slightly
- // shifted positions
- std::vector<Point<dim> > diff_points (quadrature.size());
-
- differences.resize(2*dim);
- for (unsigned int d=0; d<dim; ++d)
- {
- Tensor<1,dim> shift;
- shift[d] = fd_step_length;
-
- // generate points and FEValues
- // objects shifted in
- // plus-direction. note that
- // they only need to compute
- // gradients, not more
- for (unsigned int i=0; i<diff_points.size(); ++i)
- diff_points[i] = quadrature.point(i) + shift;
- const Quadrature<dim> plus_quad (diff_points, quadrature.get_weights());
- differences[d] = new FEValues<dim,spacedim> (mapping, *element,
- plus_quad, update_gradients);
-
- // now same in minus-direction
- for (unsigned int i=0; i<diff_points.size(); ++i)
- diff_points[i] = quadrature.point(i) - shift;
- const Quadrature<dim> minus_quad (diff_points, quadrature.get_weights());
- differences[d+dim] = new FEValues<dim,spacedim> (mapping, *element,
- minus_quad, update_gradients);
- }
-}
-
-
-
template <int dim, int spacedim>
UpdateFlags
FiniteElement<dim,spacedim>::InternalDataBase::current_update_flags () const
-template<>
-void
-FiniteElement<1,2>::compute_2nd (
- const Mapping<1,2> &,
- const Triangulation<1,2>::cell_iterator &,
- const unsigned int,
- const Mapping<1,2>::InternalDataBase &,
- const InternalDataBase &,
- internal::FEValues::FiniteElementRelatedData<1,2> &) const
-{
- Assert(false, ExcNotImplemented());
-}
-
-
-template<>
-void
-FiniteElement<1,3>::compute_2nd (
- const Mapping<1,3> &,
- const Triangulation<1,3>::cell_iterator &,
- const unsigned int,
- const Mapping<1,3>::InternalDataBase &,
- const InternalDataBase &,
- internal::FEValues::FiniteElementRelatedData<1,3> &) const
-{
- Assert(false, ExcNotImplemented());
-}
-
-
-
-template<>
-void
-FiniteElement<2,3>::compute_2nd (
- const Mapping<2,3> &,
- const Triangulation<2,3>::cell_iterator &,
- const unsigned int,
- const Mapping<2,3>::InternalDataBase &,
- const InternalDataBase &,
- internal::FEValues::FiniteElementRelatedData<2,3> &) const
-{
- Assert(false, ExcNotImplemented());
-}
-
-
-
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::compute_2nd (
- const Mapping<dim,spacedim> &mapping,
- const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- const unsigned int offset,
- const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
- const InternalDataBase &fe_internal,
- internal::FEValues::FiniteElementRelatedData<dim,spacedim> &data) const
-{
- Assert ((fe_internal.update_each | fe_internal.update_once)
- & update_hessians,
- ExcInternalError());
-
- // there is nothing to do if there are no degrees of freedom (e.g., in an
- // FE_Nothing)
- if (this->dofs_per_cell == 0)
- return;
-
-// make sure we have as many entries as there are nonzero components
-// Assert (data.shape_hessians.size() ==
-// std::accumulate (n_nonzero_components_table.begin(),
-// n_nonzero_components_table.end(),
-// 0U),
-// ExcInternalError());
- // Number of quadrature points
- const unsigned int n_q_points = data.shape_hessians[0].size();
-
- // first reinit the fe_values
- // objects used for the finite
- // differencing stuff
- for (unsigned int d=0; d<dim; ++d)
- {
- fe_internal.differences[d]->reinit(cell);
- fe_internal.differences[d+dim]->reinit(cell);
- Assert(offset <= fe_internal.differences[d]->n_quadrature_points - n_q_points,
- ExcIndexRange(offset, 0, fe_internal.differences[d]->n_quadrature_points
- - n_q_points));
- }
-
- // collection of difference
- // quotients of gradients in each
- // direction (first index) and at
- // all q-points (second index)
- std::vector<std::vector<Tensor<1,dim> > >
- diff_quot (spacedim, std::vector<Tensor<1,dim> > (n_q_points));
- std::vector<Tensor<1,spacedim> > diff_quot2 (n_q_points);
-
- // for all nonzero components of
- // all shape functions at all
- // quadrature points and difference
- // quotients in all directions:
- unsigned int total_index = 0;
- for (unsigned int shape_index=0; shape_index<this->dofs_per_cell; ++shape_index)
- for (unsigned int n=0; n<n_nonzero_components(shape_index); ++n, ++total_index)
- {
- for (unsigned int d1=0; d1<dim; ++d1)
- for (unsigned int q=0; q<n_q_points; ++q)
- {
- // get gradient at points
- // shifted slightly to
- // the right and to the
- // left in the present
- // coordinate direction
- //
- // note that things
- // might be more
- // difficult if the
- // shape function has
- // more than one
- // non-zero component,
- // so find out about
- // the actual component
- // if necessary
- Tensor<1,spacedim> right, left;
- if (is_primitive(shape_index))
- {
- right = fe_internal.differences[d1]->shape_grad(shape_index, q+offset);
- left = fe_internal.differences[d1+dim]->shape_grad(shape_index, q+offset);
- }
- else
- {
- // get the
- // component index
- // of the n-th
- // nonzero
- // component
- unsigned int component=0;
- for (unsigned int nonzero_comp=0; component<this->n_components();
- ++component)
- if (nonzero_components[shape_index][component] == true)
- {
- ++nonzero_comp;
- // check
- // whether we
- // have found
- // the
- // component
- // we are
- // looking
- // for. note
- // that
- // nonzero_comp
- // is 1-based
- // by the way
- // we compute
- // it
- if (nonzero_comp == n+1)
- break;
- }
- Assert (component < this->n_components(),
- ExcInternalError());
-
- right = fe_internal.differences[d1]
- ->shape_grad_component(shape_index, q+offset, component);
- left = fe_internal.differences[d1+dim]
- ->shape_grad_component(shape_index, q+offset, component);
- };
-
- // compute the second
- // derivative from a
- // symmetric difference
- // approximation
- for (unsigned int d=0; d<spacedim; ++d)
- diff_quot[d][q][d1] = 1./(2*fd_step_length) * (right[d]-left[d]);
- }
-
- // up to now we still have
- // difference quotients on the
- // unit cell, so transform it
- // to something on the real
- // cell
- for (unsigned int d=0; d<spacedim; ++d)
- {
- mapping.transform (diff_quot[d], mapping_covariant, mapping_internal, diff_quot2);
-
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d1=0; d1<spacedim; ++d1)
- data.shape_hessians[total_index][q][d][d1]
- = diff_quot2[q][d1];
- }
- }
-}
-
-
-
template <int dim, int spacedim>
std::vector<unsigned int>
FiniteElement<dim,spacedim>::compute_n_nonzero_components (
projects / dealii.git / commitdiff