From a325b28ab255279c95d12a616254b2830ff65de5 Mon Sep 17 00:00:00 2001 From: Maien Hamed Date: Mon, 7 Sep 2015 18:03:48 +0200 Subject: [PATCH] removed compute_2nd as hessians are analytically computed --- include/deal.II/fe/fe.h | 18 --- source/fe/fe.cc | 242 ---------------------------------------- 2 files changed, 260 deletions(-) diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index d24f96272a..ca8ee1035a 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -519,14 +519,6 @@ public: */ virtual ~InternalDataBase (); - /** - * Initialize some pointers used in the computation of second derivatives - * by finite differencing of gradients. - */ - void initialize_2nd (const FiniteElement *element, - const Mapping &mapping, - const Quadrature &quadrature); - /** * Values updated by the constructor or by reinit. */ @@ -2213,16 +2205,6 @@ protected: TableIndices<2> interface_constraints_size () const; - /** - * Compute second derivatives by finite differences of gradients. - */ - void compute_2nd (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const unsigned int offset, - const typename Mapping::InternalDataBase &mapping_internal, - const InternalDataBase &fe_internal, - dealii::internal::FEValues::FiniteElementRelatedData &data) const; - /** * Given the pattern of nonzero components for each shape function, compute * for each entry how many components are non-zero for each shape function. diff --git a/source/fe/fe.cc b/source/fe/fe.cc index 28d06a68f6..291eabc240 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -73,58 +73,6 @@ FiniteElement::InternalDataBase::memory_consumption () const - -template -void -FiniteElement:: -InternalDataBase::initialize_2nd (const FiniteElement *element, - const Mapping &mapping, - const Quadrature &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 > diff_points (quadrature.size()); - - differences.resize(2*dim); - for (unsigned int d=0; d 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 plus_quad (diff_points, quadrature.get_weights()); - differences[d] = new FEValues (mapping, *element, - plus_quad, update_gradients); - - // now same in minus-direction - for (unsigned int i=0; i minus_quad (diff_points, quadrature.get_weights()); - differences[d+dim] = new FEValues (mapping, *element, - minus_quad, update_gradients); - } -} - - - template UpdateFlags FiniteElement::InternalDataBase::current_update_flags () const @@ -1283,196 +1231,6 @@ FiniteElement::memory_consumption () 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 -void -FiniteElement::compute_2nd ( - const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const unsigned int offset, - const typename Mapping::InternalDataBase &mapping_internal, - const InternalDataBase &fe_internal, - internal::FEValues::FiniteElementRelatedData &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; dreinit(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 > > - diff_quot (spacedim, std::vector > (n_q_points)); - std::vector > 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_indexdofs_per_cell; ++shape_index) - for (unsigned int n=0; n 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; componentn_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 std::vector FiniteElement::compute_n_nonzero_components ( -- 2.39.5