*/
virtual UpdateFlags update_each (const UpdateFlags flags) const;
-
/**
* Fields of cell-independent data.
*
* multiplication) when visiting an actual cell.
*/
std::vector<std::vector<Tensor<1,dim> > > shape_gradients;
+
+ /**
+ * Array with shape function hessians in quadrature points. There is one
+ * row for each shape function, containing values for each quadrature
+ * point.
+ *
+ * We store the hessians in the quadrature points on the unit cell. We
+ * then only have to apply the transformation when visiting an actual cell.
+ */
+ std::vector<std::vector<Tensor<2,dim> > > shape_hessians;
+
+ /**
+ * Scratch array to store temporary values during hessian calculations in
+ * actual cells.
+ */
+ mutable std::vector<Tensor<2,dim> > untransformed_shape_hessians;
};
+ /**
+ * Correct the hessian in the reference cell by subtracting the term corresponding
+ * to the Jacobian gradient for one degree of freedom. The result being given by:
+ *
+ * \frac{\partial^2 \phi_i}{\partial\hat{x}_J\partial\hat{x}_K}
+ * - \frac{\partial \phi_i}{\partial {x}_l}
+ * \left( \frac{\partial^2{x}_l}{\partial\hat{x}_J\partial\hat{x}_K} \right)
+ *
+ * After this correction, the shape hessians are simply a mapping_covariant_gradient
+ * transformation.
+ */
+ void
+ correct_untransformed_hessians (VectorSlice< std::vector<Tensor<2, dim> > > uncorrected_shape_hessians,
+ const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
+ const internal::FEValues::FiniteElementRelatedData<dim,spacedim> &fevalues_data,
+ const unsigned int n_q_points,
+ const unsigned int dof) const;
+
/**
* The polynomial space. Its type is given by the template parameter POLY.
*/
if (flags & update_gradients)
out |= update_gradients | update_covariant_transformation;
if (flags & update_hessians)
- out |= update_hessians | update_covariant_transformation;
+ out |= update_hessians | update_covariant_transformation
+ | update_gradients | update_jacobian_grads;
if (flags & update_cell_normal_vectors)
out |= update_cell_normal_vectors | update_JxW_values;
template <class POLY, int dim, int spacedim>
typename FiniteElement<dim,spacedim>::InternalDataBase *
FE_Poly<POLY,dim,spacedim>::get_data (const UpdateFlags update_flags,
- const Mapping<dim,spacedim> &mapping,
+ const Mapping<dim,spacedim> &,
const Quadrature<dim> &quadrature) const
{
// generate a new data object and
// then initialize some objects for
// that
if (flags & update_hessians)
- data->initialize_2nd (this, mapping, quadrature);
+ {
+ grad_grads.resize (this->dofs_per_cell);
+ data->shape_hessians.resize (this->dofs_per_cell,
+ std::vector<Tensor<2,dim> > (n_q_points));
+ data->untransformed_shape_hessians.resize (n_q_points);
+ }
// next already fill those fields
// of which we have information by
// unit cell, and need to be
// transformed when visiting an
// actual cell
- if (flags & (update_values | update_gradients))
+ if (flags & (update_values | update_gradients | update_hessians))
for (unsigned int i=0; i<n_q_points; ++i)
{
poly_space.compute(quadrature.point(i),
if (flags & update_gradients)
for (unsigned int k=0; k<this->dofs_per_cell; ++k)
data->shape_gradients[k][i] = grads[k];
+
+ if (flags & update_hessians)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+ data->shape_hessians[k][i] = grad_grads[k];
}
return data;
}
template <class POLY, int dim, int spacedim>
void
FE_Poly<POLY,dim,spacedim>::
-fill_fe_values (const Mapping<dim,spacedim> &mapping,
- const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
+fill_fe_values (const Mapping<dim,spacedim> &mapping,
+ const typename Triangulation<dim,spacedim>::cell_iterator &,
+ const Quadrature<dim> &quadrature,
+ const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
const typename FiniteElement<dim,spacedim>::InternalDataBase &fedata,
- const internal::FEValues::MappingRelatedData<dim,spacedim> &,
- internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data,
- const CellSimilarity::Similarity cell_similarity) const
+ const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
+ internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data,
+ const CellSimilarity::Similarity cell_similarity) const
{
// convert data object to internal
// data for this class. fails with
mapping.transform(fe_data.shape_gradients[k],
output_data.shape_gradients[k],
mapping_internal, mapping_covariant);
- }
- if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
- this->compute_2nd (mapping, cell, QProjector<dim>::DataSetDescriptor::cell(),
- mapping_internal, fe_data,
- output_data);
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ {
+ // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
+ for (unsigned int i=0; i<quadrature.size(); ++i)
+ {
+ fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
+ }
+
+ correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
+ mapping_data, output_data, quadrature.size(), k);
+
+ mapping.transform(fe_data.untransformed_shape_hessians,
+ output_data.shape_hessians[k],
+ mapping_internal, mapping_covariant_gradient);
+ }
+ }
}
template <class POLY, int dim, int spacedim>
void
FE_Poly<POLY,dim,spacedim>::
-fill_fe_face_values (const Mapping<dim,spacedim> &mapping,
- const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- const unsigned int face,
- const Quadrature<dim-1> &quadrature,
+fill_fe_face_values (const Mapping<dim,spacedim> &mapping,
+ const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ const unsigned int face,
+ const Quadrature<dim-1> &quadrature,
const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
- const typename FiniteElement<dim,spacedim>::InternalDataBase &fedata,
- const internal::FEValues::MappingRelatedData<dim,spacedim> &,
- internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
+ const typename FiniteElement<dim,spacedim>::InternalDataBase &fedata,
+ const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
+ internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
{
// convert data object to internal
// data for this class. fails with
mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
output_data.shape_gradients[k],
mapping_internal, mapping_covariant);
- }
- if (flags & update_hessians)
- this->compute_2nd (mapping, cell, offset, mapping_internal, fe_data,
- output_data);
+ if (flags & update_hessians)
+ {
+ // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
+ for (unsigned int i=0; i<quadrature.size(); ++i)
+ {
+ fe_data.untransformed_shape_hessians[i+offset] = fe_data.shape_hessians[k][i+offset];
+ }
+
+ correct_untransformed_hessians(VectorSlice< std::vector<Tensor<2,dim> > >
+ ( fe_data.untransformed_shape_hessians, offset , quadrature.size()),
+ mapping_data, output_data, quadrature.size(), k);
+
+ mapping.transform(make_slice(fe_data.untransformed_shape_hessians, offset, quadrature.size()),
+ output_data.shape_hessians[k], mapping_internal, mapping_covariant_gradient);
+ }
+ }
}
-
-
template <class POLY, int dim, int spacedim>
void
FE_Poly<POLY,dim,spacedim>::
-fill_fe_subface_values (const Mapping<dim,spacedim> &mapping,
- const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- const unsigned int face,
- const unsigned int subface,
- const Quadrature<dim-1> &quadrature,
+fill_fe_subface_values (const Mapping<dim,spacedim> &mapping,
+ const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ const unsigned int face,
+ const unsigned int subface,
+ const Quadrature<dim-1> &quadrature,
const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
- const typename FiniteElement<dim,spacedim>::InternalDataBase &fedata,
- const internal::FEValues::MappingRelatedData<dim,spacedim> &,
- internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
+ const typename FiniteElement<dim,spacedim>::InternalDataBase &fedata,
+ const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
+ internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
{
// convert data object to internal
// data for this class. fails with
mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
output_data.shape_gradients[k],
mapping_internal, mapping_covariant);
- }
- if (flags & update_hessians)
- this->compute_2nd (mapping, cell, offset, mapping_internal, fe_data,
- output_data);
+ if (flags & update_hessians)
+ {
+ // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
+ for (unsigned int i=0; i<quadrature.size(); ++i)
+ {
+ fe_data.untransformed_shape_hessians[i+offset] = fe_data.shape_hessians[k][i+offset];
+ }
+
+ correct_untransformed_hessians(VectorSlice< std::vector<Tensor<2,dim> > >
+ ( fe_data.untransformed_shape_hessians, offset , quadrature.size()),
+ mapping_data, output_data, quadrature.size(), k);
+
+ mapping.transform(make_slice(fe_data.untransformed_shape_hessians, offset, quadrature.size()),
+ output_data.shape_hessians[k], mapping_internal, mapping_covariant_gradient);
+ }
+ }
}
+template <class POLY, int dim, int spacedim>
+void
+FE_Poly<POLY,dim,spacedim>::
+correct_untransformed_hessians (VectorSlice< std::vector<Tensor<2, dim> > > uncorrected_shape_hessians,
+ const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
+ const internal::FEValues::FiniteElementRelatedData<dim,spacedim> &fevalues_data,
+ const unsigned int n_q_points,
+ const unsigned int dof) const
+{
+ for (unsigned int i=0; i<n_q_points; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int l=0; l<dim; ++l)
+ for (unsigned int n=0; n<spacedim; ++n)
+ uncorrected_shape_hessians[i][j][l] -= fevalues_data.shape_gradients[dof][i][n]
+ * mapping_data.jacobian_grads[i][n][l][j];
+}
namespace internal
{
public:
/**
* Constructor. Is called by the @p get_data function. Sets the size of
- * the @p base_fe_datas vector to @p n_base_elements and initializes the
- * compute_hessians field.
+ * the @p base_fe_datas vector to @p n_base_elements.
*/
- InternalData (const unsigned int n_base_elements,
- const bool compute_hessians);
+ InternalData (const unsigned int n_base_elements);
/**
* Destructor. Deletes all @p InternalDatas whose pointers are stored by
*/
~InternalData();
- /**
- * Flag indicating whether second derivatives shall be computed.
- */
- const bool compute_hessians;
-
/**
* Gives write-access to the pointer to a @p InternalData of the @p
* base_noth base element.
#include <deal.II/base/polynomials_piecewise.h>
#include <deal.II/fe/fe_poly.h>
#include <deal.II/fe/fe_values.h>
-
#include <deal.II/fe/fe_poly.templates.h>
DEAL_II_NAMESPACE_OPEN
template <>
void
FE_Poly<TensorProductPolynomials<1>,1,2>::
-fill_fe_values (const Mapping<1,2> &mapping,
- const Triangulation<1,2>::cell_iterator &cell,
- const Quadrature<1> &quadrature,
- const Mapping<1,2>::InternalDataBase &mapping_internal,
- const FiniteElement<1,2>::InternalDataBase &fedata,
- const internal::FEValues::MappingRelatedData<1,2> &,
+fill_fe_values (const Mapping<1,2> &mapping,
+ const Triangulation<1,2>::cell_iterator &,
+ const Quadrature<1> &quadrature,
+ const Mapping<1,2>::InternalDataBase &mapping_internal,
+ const FiniteElement<1,2>::InternalDataBase &fedata,
+ const internal::FEValues::MappingRelatedData<1,2> &mapping_data,
internal::FEValues::FiniteElementRelatedData<1,2> &output_data,
- const CellSimilarity::Similarity cell_similarity) const
+ const CellSimilarity::Similarity cell_similarity) const
{
// convert data object to internal
// data for this class. fails with
mapping.transform(fe_data.shape_gradients[k],
output_data.shape_gradients[k],
mapping_internal, mapping_covariant);
- }
- if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
- this->compute_2nd (mapping, cell, QProjector<1>::DataSetDescriptor::cell(),
- mapping_internal, fe_data,
- output_data);
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ {
+ // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
+ for (unsigned int i=0; i<quadrature.size(); ++i)
+ {
+ fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
+ }
+
+ correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
+ mapping_data, output_data, quadrature.size(), k);
+
+ mapping.transform(fe_data.untransformed_shape_hessians,
+ output_data.shape_hessians[k],
+ mapping_internal, mapping_covariant_gradient);
+ }
+ }
}
template <>
void
FE_Poly<TensorProductPolynomials<2>,2,3>::
-fill_fe_values (const Mapping<2,3> &mapping,
- const Triangulation<2,3>::cell_iterator &cell,
- const Quadrature<2> &quadrature,
- const Mapping<2,3>::InternalDataBase &mapping_internal,
- const FiniteElement<2,3>::InternalDataBase &fedata,
- const internal::FEValues::MappingRelatedData<2,3> &,
+fill_fe_values (const Mapping<2,3> &mapping,
+ const Triangulation<2,3>::cell_iterator &,
+ const Quadrature<2> &quadrature,
+ const Mapping<2,3>::InternalDataBase &mapping_internal,
+ const FiniteElement<2,3>::InternalDataBase &fedata,
+ const internal::FEValues::MappingRelatedData<2,3> &mapping_data,
internal::FEValues::FiniteElementRelatedData<2,3> &output_data,
- const CellSimilarity::Similarity cell_similarity) const
+ const CellSimilarity::Similarity cell_similarity) const
{
// assert that the following dynamics
mapping.transform(fe_data.shape_gradients[k],
output_data.shape_gradients[k],
mapping_internal, mapping_covariant);
- }
- if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
- this->compute_2nd (mapping, cell, QProjector<2>::DataSetDescriptor::cell(),
- mapping_internal, fe_data,
- output_data);
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ {
+ // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
+ for (unsigned int i=0; i<quadrature.size(); ++i)
+ {
+ fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
+ }
+
+ correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
+ mapping_data, output_data, quadrature.size(), k);
+
+ mapping.transform(fe_data.untransformed_shape_hessians,
+ output_data.shape_hessians[k],
+ mapping_internal, mapping_covariant_gradient);
+ }
+ }
}
template <>
void
FE_Poly<PolynomialSpace<1>,1,2>::
-fill_fe_values (const Mapping<1,2> &mapping,
- const Triangulation<1,2>::cell_iterator &cell,
- const Quadrature<1> &quadrature,
- const Mapping<1,2>::InternalDataBase &mapping_internal,
- const FiniteElement<1,2>::InternalDataBase &fedata,
- const internal::FEValues::MappingRelatedData<1,2> &,
+fill_fe_values (const Mapping<1,2> &mapping,
+ const Triangulation<1,2>::cell_iterator &,
+ const Quadrature<1> &quadrature,
+ const Mapping<1,2>::InternalDataBase &mapping_internal,
+ const FiniteElement<1,2>::InternalDataBase &fedata,
+ const internal::FEValues::MappingRelatedData<1,2> &mapping_data,
internal::FEValues::FiniteElementRelatedData<1,2> &output_data,
- const CellSimilarity::Similarity cell_similarity) const
+ const CellSimilarity::Similarity cell_similarity) const
{
// convert data object to internal
// data for this class. fails with
mapping.transform(fe_data.shape_gradients[k],
output_data.shape_gradients[k],
mapping_internal, mapping_covariant);
- }
- if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
- this->compute_2nd (mapping, cell, QProjector<1>::DataSetDescriptor::cell(),
- mapping_internal, fe_data,
- output_data);
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ {
+ // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
+ for (unsigned int i=0; i<quadrature.size(); ++i)
+ {
+ fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
+ }
+
+ correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
+ mapping_data, output_data, quadrature.size(), k);
+
+ mapping.transform(fe_data.untransformed_shape_hessians,
+ output_data.shape_hessians[k],
+ mapping_internal, mapping_covariant_gradient);
+ }
+ }
}
template <>
void
FE_Poly<PolynomialSpace<2>,2,3>::
-fill_fe_values (const Mapping<2,3> &mapping,
- const Triangulation<2,3>::cell_iterator &cell,
- const Quadrature<2> &quadrature,
- const Mapping<2,3>::InternalDataBase &mapping_internal,
- const FiniteElement<2,3>::InternalDataBase &fedata,
- const internal::FEValues::MappingRelatedData<2,3> &,
+fill_fe_values (const Mapping<2,3> &mapping,
+ const Triangulation<2,3>::cell_iterator &,
+ const Quadrature<2> &quadrature,
+ const Mapping<2,3>::InternalDataBase &mapping_internal,
+ const FiniteElement<2,3>::InternalDataBase &fedata,
+ const internal::FEValues::MappingRelatedData<2,3> &mapping_data,
internal::FEValues::FiniteElementRelatedData<2,3> &output_data,
- const CellSimilarity::Similarity cell_similarity) const
+ const CellSimilarity::Similarity cell_similarity) const
{
Assert (dynamic_cast<const InternalData *> (&fedata) != 0, ExcInternalError());
const InternalData &fe_data = static_cast<const InternalData &> (fedata);
mapping.transform(fe_data.shape_gradients[k],
output_data.shape_gradients[k],
mapping_internal, mapping_covariant);
- }
- if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
- this->compute_2nd (mapping, cell, QProjector<2>::DataSetDescriptor::cell(),
- mapping_internal, fe_data,
- output_data);
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ {
+ // compute the hessians in the unit cell (accounting for the Jacobian gradiant)
+ for (unsigned int i=0; i<quadrature.size(); ++i)
+ {
+ fe_data.untransformed_shape_hessians[i] = fe_data.shape_hessians[k][i];
+ }
+
+ correct_untransformed_hessians (fe_data.untransformed_shape_hessians,
+ mapping_data, output_data, quadrature.size(), k);
+
+ mapping.transform(fe_data.untransformed_shape_hessians,
+ output_data.shape_hessians[k],
+ mapping_internal, mapping_covariant_gradient);
+ }
+ }
}
template <int dim, int spacedim>
-FESystem<dim,spacedim>::InternalData::InternalData(const unsigned int n_base_elements,
- const bool compute_hessians)
+FESystem<dim,spacedim>::InternalData::InternalData(const unsigned int n_base_elements)
:
- compute_hessians (compute_hessians),
base_fe_datas(n_base_elements),
base_fe_output_objects(n_base_elements)
{}
for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
out |= base_element(base_no).update_each(flags);
- // second derivatives are handled
- // by the top-level finite element,
- // rather than by the base elements
- // since it is generated by finite
- // differencing. if second
- // derivatives are requested, we
- // therefore have to set the
- // respective flag since the base
- // elements don't have them
- if (flags & update_hessians)
- out |= update_hessians | update_covariant_transformation;
-
return out;
}
const Quadrature<dim> &quadrature) const
{
UpdateFlags flags = flags_;
- InternalData *data = new InternalData(this->n_base_elements(),
- flags & update_hessians);
+ InternalData *data = new InternalData(this->n_base_elements());
data->update_once = update_once (flags);
data->update_each = update_each (flags);
flags = data->update_once | data->update_each;
- UpdateFlags sub_flags = flags;
- // if second derivatives through
- // finite differencing are required,
- // then initialize some objects for
- // that
- if (data->compute_hessians)
- {
- // delete
- // update_hessians
- // from flags list
- sub_flags = UpdateFlags (sub_flags ^ update_hessians);
- data->initialize_2nd (this, mapping, quadrature);
- }
-
// get data objects from each of
// the base elements and store them
for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
// create internal objects for base elements and initialize their
// respective output objects
typename FiniteElement<dim,spacedim>::InternalDataBase *base_fe_data =
- base_element(base_no).get_data(sub_flags, mapping, quadrature);
+ base_element(base_no).get_data(flags, mapping, quadrature);
internal::FEValues::FiniteElementRelatedData<dim,spacedim> &base_fe_output_object
= data->get_fe_output_object(base_no);
// then store the pointer to the base internal object
data->set_fe_data(base_no, base_fe_data);
-
- // make sure that *we* compute
- // second derivatives, base
- // elements should not do it
- Assert (!(base_fe_data->update_each & update_hessians),
- ExcInternalError());
- Assert (!(base_fe_data->update_once & update_hessians),
- ExcInternalError());
}
data->update_flags = data->update_once |
data->update_each;
const Quadrature<dim-1> &quadrature) const
{
UpdateFlags flags = flags_;
- InternalData *data = new InternalData(this->n_base_elements(),
- flags & update_hessians);
+ InternalData *data = new InternalData(this->n_base_elements());
data->update_once = update_once (flags);
data->update_each = update_each (flags);
flags = data->update_once | data->update_each;
- UpdateFlags sub_flags = flags;
- if (data->compute_hessians)
- {
- sub_flags = UpdateFlags (sub_flags ^ update_hessians);
- data->initialize_2nd (this, mapping, QProjector<dim>::project_to_all_faces(quadrature));
- }
-
for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
{
// create internal objects for base elements and initialize their
// respective output objects
typename FiniteElement<dim,spacedim>::InternalDataBase *base_fe_data =
- base_element(base_no).get_face_data(sub_flags, mapping, quadrature);
+ base_element(base_no).get_face_data(flags, mapping, quadrature);
internal::FEValues::FiniteElementRelatedData<dim,spacedim> &base_fe_output_object
= data->get_fe_output_object(base_no);
// then store the pointer to the base internal object
data->set_fe_data(base_no, base_fe_data);
-
- Assert (!(base_fe_data->update_each & update_hessians),
- ExcInternalError());
- Assert (!(base_fe_data->update_once & update_hessians),
- ExcInternalError());
}
data->update_flags = data->update_once |
data->update_each;
const Quadrature<dim-1> &quadrature) const
{
UpdateFlags flags = flags_;
- InternalData *data = new InternalData(this->n_base_elements(),
- flags & update_hessians);
+ InternalData *data = new InternalData(this->n_base_elements());
data->update_once = update_once (flags);
data->update_each = update_each (flags);
flags = data->update_once | data->update_each;
- UpdateFlags sub_flags = flags;
- if (data->compute_hessians)
- {
- sub_flags = UpdateFlags (sub_flags ^ update_hessians);
- data->initialize_2nd (this, mapping, QProjector<dim>::project_to_all_subfaces(quadrature));
- }
-
for (unsigned int base_no=0; base_no<this->n_base_elements(); ++base_no)
{
// create internal objects for base elements and initialize their
// respective output objects
typename FiniteElement<dim,spacedim>::InternalDataBase *base_fe_data =
- base_element(base_no).get_subface_data(sub_flags, mapping, quadrature);
+ base_element(base_no).get_subface_data(flags, mapping, quadrature);
internal::FEValues::FiniteElementRelatedData<dim,spacedim> &base_fe_output_object
= data->get_fe_output_object(base_no);
// then store the pointer to the base internal object
data->set_fe_data(base_no, base_fe_data);
-
- Assert (!(base_fe_data->update_each & update_hessians),
- ExcInternalError());
- Assert (!(base_fe_data->update_once & update_hessians),
- ExcInternalError());
}
data->update_flags = data->update_once |
data->update_each;
output_data.shape_gradients[out_index+s][q] =
base_data.shape_gradients[in_index+s][q];
- // _we_ handle computation of second derivatives, so the
- // base elements should not have computed them!
- Assert (!(base_flags & update_hessians),
- ExcInternalError());
+ if (base_flags & update_hessians)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ output_data.shape_hessians[out_index+s][q] =
+ base_data.shape_hessians[in_index+s][q];
+
}
}
}
const unsigned int face_no,
const unsigned int sub_no,
const Quadrature<dim_1> &quadrature,
- const CellSimilarity::Similarity cell_similarity,
+ const CellSimilarity::Similarity ,
const typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
const typename FiniteElement<dim,spacedim>::InternalDataBase &fedata,
const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,
internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const
{
- const unsigned int n_q_points = quadrature.size();
-
// convert data object to internal
// data for this class. fails with
// an exception if that is not
fe_data.update_flags);
- if (flags & (update_values | update_gradients))
+ if (flags & (update_values | update_gradients | update_hessians))
{
// let base elements update the necessary data
Threads::TaskGroup<> task_group;
std_cxx11::ref(output_data))));
task_group.join_all();
}
-
- if (fe_data.compute_hessians && cell_similarity != CellSimilarity::translation)
- {
- unsigned int offset = 0;
- if (face_no != invalid_face_number)
- {
- if (sub_no == invalid_face_number)
- offset=QProjector<dim>::DataSetDescriptor
- ::face(face_no,
- cell->face_orientation(face_no),
- cell->face_flip(face_no),
- cell->face_rotation(face_no),
- n_q_points);
- else
- offset=QProjector<dim>::DataSetDescriptor
- ::subface(face_no, sub_no,
- cell->face_orientation(face_no),
- cell->face_flip(face_no),
- cell->face_rotation(face_no),
- n_q_points,
- cell->subface_case(face_no));
- }
-
- this->compute_2nd (mapping, cell, offset, mapping_internal, fe_data,
- output_data);
- }
}
}
case mapping_piola:
+ {
+ Assert (data.update_flags & update_contravariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+ Assert (data.update_flags & update_volume_elements,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
+
+ for (unsigned int i=0; i<output.size(); ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ output[i][d1][d2] = input[i][d1][d2] * data.length[d2]
+ / data.volume_element;
+ return;
+ }
+
+ case mapping_piola_gradient:
{
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
switch (mapping_type)
{
+ case mapping_covariant:
+ {
+ Assert (data.update_flags & update_covariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+ for (unsigned int i=0; i<output.size(); ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ output[i][d1][d2] = input[i][d1][d2] / data.length[d2];
+ return;
+ }
+
+ case mapping_contravariant:
+ {
+ Assert (data.update_flags & update_contravariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+
+ for (unsigned int i=0; i<output.size(); ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ output[i][d1][d2] = input[i][d1][d2] * data.length[d2];
+ return;
+ }
+
+ case mapping_covariant_gradient:
+ {
+ Assert (data.update_flags & update_covariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_covariant_transformation"));
+
+ for (unsigned int i=0; i<output.size(); ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ output[i][d1][d2] = input[i][d1][d2] / data.length[d2] / data.length[d1];
+ return;
+ }
+
+ case mapping_contravariant_gradient:
+ {
+ Assert (data.update_flags & update_contravariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+
+ for (unsigned int i=0; i<output.size(); ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ output[i][d1][d2] = input[i][d1][d2] * data.length[d2] / data.length[d1];
+ return;
+ }
+
+ case mapping_piola:
+ {
+ Assert (data.update_flags & update_contravariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+ Assert (data.update_flags & update_volume_elements,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
+
+ for (unsigned int i=0; i<output.size(); ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ output[i][d1][d2] = input[i][d1][d2] * data.length[d2]
+ / data.volume_element;
+ return;
+ }
case mapping_piola_gradient:
{
Assert (data.update_flags & update_contravariant_transformation,
projects / dealii.git / commitdiff