From 5d90c842817dc0f4775e3a4500a8e94f9822f634 Mon Sep 17 00:00:00 2001 From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Thu, 8 Jan 2009 00:00:16 +0000 Subject: [PATCH] transformations moved into Mapping git-svn-id: https://svn.dealii.org/trunk@18127 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_poly.templates.h | 67 +-- .../deal.II/include/fe/fe_raviart_thomas.h | 4 +- deal.II/deal.II/include/fe/mapping.h | 50 ++- .../deal.II/include/fe/mapping_cartesian.h | 27 +- deal.II/deal.II/include/fe/mapping_q.h | 40 -- deal.II/deal.II/include/fe/mapping_q1.h | 24 -- deal.II/deal.II/source/fe/fe_poly_tensor.cc | 283 +++++-------- .../deal.II/source/fe/fe_raviart_thomas.cc | 20 +- .../source/fe/fe_raviart_thomas_nodal.cc | 29 +- deal.II/deal.II/source/fe/mapping.cc | 65 ++- .../deal.II/source/fe/mapping_cartesian.cc | 252 +++++------- deal.II/deal.II/source/fe/mapping_q.cc | 212 ++-------- deal.II/deal.II/source/fe/mapping_q1.cc | 381 +++++++----------- 13 files changed, 546 insertions(+), 908 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_poly.templates.h b/deal.II/deal.II/include/fe/fe_poly.templates.h index b579d00b4e..337cf7fbf1 100644 --- a/deal.II/deal.II/include/fe/fe_poly.templates.h +++ b/deal.II/deal.II/include/fe/fe_poly.templates.h @@ -260,14 +260,8 @@ FE_Poly<TensorProductPolynomials<1>,1,2>::fill_fe_values (const Mapping<1,2> &ma data.shape_values(k,i) = fe_data.shape_values[k][i]; if (flags & update_gradients) - { - Assert (data.shape_gradients[k].size() <= - fe_data.shape_gradients[k].size(), - ExcInternalError()); - mapping.transform(fe_data.shape_gradients[k], - data.shape_gradients[k], - mapping_data, mapping_covariant); - } + mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k], + mapping_data, mapping_covariant); } // const typename QProjector<1>::DataSetDescriptor dsd; @@ -303,14 +297,8 @@ FE_Poly<TensorProductPolynomials<2>,2,3>::fill_fe_values (const Mapping<2,3> &ma data.shape_values(k,i) = fe_data.shape_values[k][i]; if (flags & update_gradients) - { - Assert (data.shape_gradients[k].size() <= - fe_data.shape_gradients[k].size(), - ExcInternalError()); - mapping.transform(fe_data.shape_gradients[k], - data.shape_gradients[k], + mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k], mapping_data, mapping_covariant); - } } /* const QProjector<2>::DataSetDescriptor dsd; */ @@ -327,12 +315,13 @@ FE_Poly<TensorProductPolynomials<2>,2,3>::fill_fe_values (const Mapping<2,3> &ma template <> inline 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, - Mapping<1,2>::InternalDataBase &/* mapping_data */, - Mapping<1,2>::InternalDataBase &fedata, - FEValuesData<1,2> &data) const +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, + Mapping<1,2>::InternalDataBase& mapping_data, + Mapping<1,2>::InternalDataBase& fedata, + FEValuesData<1,2>& data) const { // convert data object to internal // data for this class. fails with @@ -353,12 +342,8 @@ FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_values (const Mapping<1,2> &/* mapping // TODO: I would think this should work. Guido if (flags & update_gradients) - { - AssertThrow(false, ExcNotImplemented()); -// mapping.transform(fe_data.shape_gradients[k], 0, -// data.shape_gradients[k], -// mapping_data, mapping_covariant); - } + mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k], + mapping_data, mapping_covariant); } // const typename QProjector<1>::DataSetDescriptor dsd; @@ -374,10 +359,10 @@ FE_Poly<PolynomialSpace<1>,1,2>::fill_fe_values (const Mapping<1,2> &/* mapping template <> inline void -FE_Poly<PolynomialSpace<2>,2,3>::fill_fe_values (const Mapping<2,3> &/* mapping */, +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, - Mapping<2,3>::InternalDataBase &/* mapping_data */, + const Quadrature<2>& quadrature, + Mapping<2,3>::InternalDataBase& mapping_data, Mapping<2,3>::InternalDataBase &fedata, FEValuesData<2,3> &data) const { @@ -394,15 +379,8 @@ FE_Poly<PolynomialSpace<2>,2,3>::fill_fe_values (const Mapping<2,3> &/* mapping if (flags & update_gradients) - { - AssertThrow(false, ExcNotImplemented()); -/* Assert (data.shape_gradients[k].size() */ -/* fe_data.shape_gradients[k].size(), */ -/* ExcInternalError()); */ -/* mapping.transform(fe_data.shape_gradients[k], 0, */ -/* data.shape_gradients[k], */ -/* mapping_data, mapping_covariant); */ - } + mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k], + mapping_data, mapping_covariant); } /* const QProjector<2>::DataSetDescriptor dsd; */ @@ -670,14 +648,9 @@ FE_Poly<POLY,dim,spacedim>::fill_fe_subface_values (const Mapping<dim,spacedim> data.shape_values(k,i) = fe_data.shape_values[k][i+offset]; if (flags & update_gradients) - { - Assert (data.shape_gradients[k].size() + offset <= - fe_data.shape_gradients[k].size(), - ExcInternalError()); - mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()), - data.shape_gradients[k], - mapping_data, mapping_covariant); - } + mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()), + data.shape_gradients[k], + mapping_data, mapping_covariant); } if (flags & update_hessians) diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h index 6895cfc790..fdb283259d 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -442,6 +442,8 @@ class FE_RaviartThomasNodal virtual FiniteElementDomination::Domination compare_for_face_domination (const FiniteElement<dim> &fe_other) const; + virtual UpdateFlags update_once (const UpdateFlags flags) const; + virtual UpdateFlags update_each (const UpdateFlags flags) const; private: /** * Only for internal use. Its diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index 51023ab054..095d4e6204 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -64,10 +64,43 @@ enum MappingType mapping_covariant = 0x0001, /// Contravariant mapping mapping_contravariant = 0x0002, +/// Mapping of the gradient of a covariant vector field + mapping_covariant_gradient = 0x0003, +/// Mapping of the gradient of a contravariant vector field + mapping_contravariant_gradient = 0x0004, /// The Piola transform usually used for Hdiv elements + /** + * Piola transform is the + * standard transformation of + * vector valued elements in + * H<sup>div</sup>. It amounts + * to a contravariant + * transformation scaled by the + * inverse of the volume + * element. + * + * If applied to a rank 2 + * tensor, the mapping class + * will apply the correct + * transformation for the + * gradient of such a vector + * field. + */ mapping_piola = 0x0100, -/// The Piola transform for the derivative of an Hdiv element - mapping_piola_gradient = 0x0101, +/// The mapping used for Nedelec elements + /** + * curl-conforming elements are + * mapped as covariant + * vectors. Nevertheless, we + * introduce a separate mapping + * type, such that we can use + * the same flag for the vector + * and its gradient (the + * transformation of the + * gradient differs from the + * one used by #mapping_covariant). + */ + mapping_nedelec = 0x0200, /// The mapping used for Raviart-Thomas elements mapping_raviart_thomas = mapping_piola, /// The mapping used for BDM elements @@ -332,42 +365,39 @@ class Mapping : public Subscriptor /** * @deprecated Use transform() instead. */ - virtual void transform_covariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input, const unsigned int offset, VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const InternalDataBase &internal) const = 0; + const InternalDataBase &internal) const; /** * @deprecated Use transform() instead. */ - virtual void transform_covariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input, const unsigned int offset, VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const InternalDataBase &internal) const = 0; + const InternalDataBase &internal) const; /** * @deprecated Use transform() instead. */ - virtual void transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input, const unsigned int offset, VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const = 0; + const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; /** * @deprecated Use transform() instead. */ - virtual void + void transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > intput, const unsigned int offset, const VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const = 0; + const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; /** * The transformed (generalized) diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index 3056a686c7..e15022238b 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -102,33 +102,8 @@ class MappingCartesian : public Mapping<dim,spacedim> transform (const VectorSlice<const std::vector<Tensor<2,dim> > > input, VectorSlice<std::vector<Tensor<2,spacedim> > > output, const typename Mapping<dim,spacedim>::InternalDataBase &internal, - const MappingType type) const; - - - virtual void - transform_covariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &internal) const; - - virtual void - transform_covariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &internal) const; - - virtual void - transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &internal) const; + const MappingType type) const; - virtual void - transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &internal) const; - virtual Point<spacedim> transform_unit_to_real_cell ( const typename Triangulation<dim,spacedim>::cell_iterator &cell, diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 288f12ce32..13e7a4546a 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -128,46 +128,6 @@ class MappingQ : public MappingQ1<dim,spacedim> const typename Mapping<dim,spacedim>::InternalDataBase &internal, const MappingType type) const; - /** - * Implementation of the - * interface in Mapping. - */ - virtual void - transform_covariant (const VectorSlice<const std::vector<Tensor<1,spacedim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - - /** - * Implementation of the - * interface in Mapping. - */ - virtual void - transform_covariant (const VectorSlice<const std::vector<Tensor<2,spacedim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - - /** - * Implementation of the - * interface in Mapping. - */ - virtual void - transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - - /** - * Implementation of the - * interface in Mapping. - */ - virtual void - transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - /** * Return the degree of the * mapping, i.e. the value which diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index fd9a81d643..0e059bb1dc 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -89,30 +89,6 @@ class MappingQ1 : public Mapping<dim,spacedim> const typename Mapping<dim,spacedim>::InternalDataBase &internal, const MappingType type) const; - virtual void - transform_covariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - - virtual void - transform_covariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - - virtual void - transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - - virtual void - transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal) const; - /** * Return a pointer to a copy of the * present object. The caller of this diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index 3ca14a3a4e..1ce943a7c7 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -407,47 +407,27 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( case mapping_covariant: case mapping_contravariant: { - // Use auxiliary vector for - // transformation std::vector<Tensor<1,dim> > shape_values (n_q_points); - if (mapping_type == mapping_covariant) - mapping.transform(fe_data.shape_values[i], - shape_values, mapping_data, mapping_covariant); - else - mapping.transform(fe_data.shape_values[i], - shape_values, mapping_data, mapping_contravariant); - - // then copy over to target: + mapping.transform(fe_data.shape_values[i], + shape_values, mapping_data, mapping_type); + for (unsigned int k=0; k<n_q_points; ++k) - { -//TODO: [GK] Remove this after change to Piola complete - - // Recompute determinant - - const double - J = (mapping_type == mapping_contravariant ? - data.JxW_values[k] / quadrature.weight(k) - : - 1.0); - - for (unsigned int d=0; d<dim; ++d) - data.shape_values(first+d,k) - = sign_change[i] * (shape_values[k][d] / J); - } + for (unsigned int d=0; d<dim; ++d) + data.shape_values(first+d,k) = shape_values[k][d]; + break; } case mapping_piola: - if (true) - { - std::vector<Tensor<1,dim> > shape_values (n_q_points); - mapping.transform(fe_data.shape_values[i], shape_values, - mapping_data, mapping_piola); - for (unsigned int k=0; k<n_q_points; ++k) - for (unsigned int d=0; d<dim; ++d) - data.shape_values(first+d,k) - = sign_change[i] * shape_values[k][d]; - break; - } + { + std::vector<Tensor<1,dim> > shape_values (n_q_points); + mapping.transform(fe_data.shape_values[i], shape_values, + mapping_data, mapping_piola); + for (unsigned int k=0; k<n_q_points; ++k) + for (unsigned int d=0; d<dim; ++d) + data.shape_values(first+d,k) + = sign_change[i] * shape_values[k][d]; + break; + } default: Assert(false, ExcNotImplemented()); @@ -470,47 +450,40 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( data.shape_gradients[first+d][k] = shape_grads1[k][d]; break; } - case mapping_covariant: { mapping.transform(fe_data.shape_grads[i], shape_grads1, - mapping_data, mapping_covariant); - for (unsigned int q=0; q<n_q_points; ++q) - shape_grads2[q] = transpose(shape_grads1[q]); - // do second transformation - mapping.transform(shape_grads2, shape_grads1, - mapping_data, mapping_covariant); - // transpose back - for (unsigned int q=0; q<n_q_points; ++q) - shape_grads2[q] = transpose(shape_grads1[q]); + mapping_data, mapping_covariant_gradient); + + for (unsigned int k=0; k<n_q_points; ++k) + for (unsigned int d=0; d<dim; ++d) + data.shape_gradients[first+d][k] = shape_grads1[k][d]; + + break; + } + case mapping_contravariant: + { + mapping.transform(fe_data.shape_grads[i], shape_grads1, + mapping_data, mapping_contravariant_gradient); for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) data.shape_gradients[first+d][k] = shape_grads2[k][d]; - + break; } case mapping_piola: - case mapping_contravariant: { mapping.transform(fe_data.shape_grads[i], shape_grads1, - mapping_data, mapping_covariant); + mapping_data, mapping_piola); - mapping.transform(shape_grads1, shape_grads2, - mapping_data, mapping_contravariant); - for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) - { - // Recompute determinant - const double J = data.JxW_values[k] / quadrature.weight(k); - data.shape_gradients[first+d][k] = sign_change[i] * - shape_grads2[k][d] / J; - } + data.shape_gradients[first+d][k] = sign_change[i] * shape_grads1[k][d]; break; } - + default: Assert(false, ExcNotImplemented()); } @@ -593,52 +566,26 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( // Use auxiliary vector // for transformation std::vector<Tensor<1,dim> > shape_values (n_q_points); - if (mapping_type == mapping_covariant) - mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_data, mapping_covariant); - else - mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_data, mapping_contravariant); + mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), + shape_values, mapping_data, mapping_type); - // then copy over to target: for (unsigned int k=0; k<n_q_points; ++k) - { - // recompute - // determinant. note that - // we have to take the - // *cell* transformation - // (cell_JxW_values) - // because the usual - // JxW_values field - // contains the - // transformation of the - // face in - // fill_fe_face_values - const double - J = (mapping_type == mapping_contravariant ? - data.cell_JxW_values[k] / quadrature.weight(k) - : - 1.0); - - for (unsigned int d=0; d<dim; ++d) - data.shape_values(first+d,k) - = sign_change[i] * (shape_values[k][d] / J); - } + for (unsigned int d=0; d<dim; ++d) + data.shape_values(first+d,k) = shape_values[k][d]; break; } case mapping_piola: - if (true) - { - std::vector<Tensor<1,dim> > shape_values (n_q_points); - mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_data, mapping_piola); - for (unsigned int k=0; k<n_q_points; ++k) - for (unsigned int d=0; d<dim; ++d) - data.shape_values(first+d,k) - = sign_change[i] * shape_values[k][d]; - break; - } + { + std::vector<Tensor<1,dim> > shape_values (n_q_points); + mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), + shape_values, mapping_data, mapping_piola); + for (unsigned int k=0; k<n_q_points; ++k) + for (unsigned int d=0; d<dim; ++d) + data.shape_values(first+d,k) + = sign_change[i] * shape_values[k][d]; + break; + } default: Assert(false, ExcNotImplemented()); @@ -665,40 +612,33 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( case mapping_covariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), - shape_grads1, mapping_data, mapping_covariant); - for (unsigned int q=0; q<n_q_points; ++q) - shape_grads2[q] = transpose(shape_grads1[q]); - // do second transformation - mapping.transform(shape_grads2, shape_grads1, - mapping_data, mapping_covariant); - // transpose back - for (unsigned int q=0; q<n_q_points; ++q) - shape_grads2[q] = transpose(shape_grads1[q]); + shape_grads1, mapping_data, mapping_covariant_gradient); for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) data.shape_gradients[first+d][k] = shape_grads1[k][d]; break; } - case mapping_piola: case mapping_contravariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), - shape_grads1, mapping_data, mapping_covariant); - - mapping.transform(shape_grads1, shape_grads2, - mapping_data, mapping_contravariant); + shape_grads1, mapping_data, mapping_contravariant_gradient); for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) - { - // recompute - // determinant in the - // same way as above - const double J = data.cell_JxW_values[k] / quadrature.weight(k); - data.shape_gradients[first+d][k] - = sign_change[i] * shape_grads2[k][d] / J; - } + data.shape_gradients[first+d][k] = shape_grads1[k][d]; + + break; + } + case mapping_piola: + { + mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), + shape_grads1, mapping_data, mapping_piola); + + for (unsigned int k=0; k<n_q_points; ++k) + for (unsigned int d=0; d<dim; ++d) + data.shape_gradients[first+d][k] = sign_change[i] * shape_grads1[k][d]; + break; } @@ -783,52 +723,26 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( // Use auxiliary vector for // transformation std::vector<Tensor<1,dim> > shape_values (n_q_points); - if (mapping_type == mapping_covariant) - mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_data, mapping_covariant); - else - mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_data, mapping_contravariant); - - // then copy over to target: + mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), + shape_values, mapping_data, mapping_type); + for (unsigned int k=0; k<n_q_points; ++k) - { - // recompute - // determinant. note that - // we have to take the - // *cell* transformation - // (cell_JxW_values) - // because the usual - // JxW_values field - // contains the - // transformation of the - // face in - // fill_fe_face_values - const double - J = (mapping_type == mapping_contravariant ? - data.cell_JxW_values[k] / quadrature.weight(k) - : - 1.0); - for (unsigned int d=0; d<dim; ++d) - data.shape_values(first+d,k) - = sign_change[i] * (shape_values[k][d] / J); - } - + data.shape_values(first+d,k) = shape_values[k][d]; + break; } case mapping_piola: - if (true) - { - std::vector<Tensor<1,dim> > shape_values (n_q_points); - mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_data, mapping_piola); - for (unsigned int k=0; k<n_q_points; ++k) - for (unsigned int d=0; d<dim; ++d) - data.shape_values(first+d,k) - = sign_change[i] * shape_values[k][d]; - break; - } + { + std::vector<Tensor<1,dim> > shape_values (n_q_points); + mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), + shape_values, mapping_data, mapping_piola); + for (unsigned int k=0; k<n_q_points; ++k) + for (unsigned int d=0; d<dim; ++d) + data.shape_values(first+d,k) + = sign_change[i] * shape_values[k][d]; + break; + } default: Assert(false, ExcNotImplemented()); } @@ -854,46 +768,41 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( case mapping_covariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), - shape_grads1, mapping_data, mapping_covariant); - for (unsigned int q=0; q<n_q_points; ++q) - shape_grads2[q] = transpose(shape_grads1[q]); - // do second transformation - mapping.transform(shape_grads2, shape_grads1, - mapping_data, mapping_covariant); - // transpose back - for (unsigned int q=0; q<n_q_points; ++q) - shape_grads2[q] = transpose(shape_grads1[q]); - + shape_grads1, mapping_data, mapping_covariant_gradient); + for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) data.shape_gradients[first+d][k] = shape_grads1[k][d]; break; } - case mapping_piola: case mapping_contravariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), - shape_grads1, mapping_data, mapping_covariant); - - mapping.transform(shape_grads1, shape_grads2, - mapping_data, mapping_contravariant); - + shape_grads1, mapping_data, mapping_contravariant_gradient); + for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) - { - // recompute - // determinant in the - // same way as above - const double J = data.cell_JxW_values[k] / quadrature.weight(k); - data.shape_gradients[first+d][k] - = sign_change[i] * shape_grads2[k][d] / J; - } + data.shape_gradients[first+d][k] = shape_grads1[k][d]; + break; } + case mapping_piola: + { + mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), + shape_grads1, + mapping_data, mapping_piola); + + for (unsigned int k=0; k<n_q_points; ++k) + for (unsigned int d=0; d<dim; ++d) + data.shape_gradients[first+d][k] = sign_change[i] * shape_grads1[k][d]; + + break; + } + default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } } } diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index 900ebb1c15..cce73d0fcc 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -50,7 +50,7 @@ FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg) Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_piola; + this->mapping_type = mapping_raviart_thomas; // First, initialize the // generalized support points and // quadrature weights, since they @@ -463,20 +463,14 @@ FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const UpdateFlags out = update_default; if (flags & update_values) - out |= update_values | update_covariant_transformation - | update_piola - | update_cell_JxW_values - | update_JxW_values; + out |= update_values | update_piola; + if (flags & update_gradients) - out |= update_gradients | update_covariant_transformation - | update_piola - | update_cell_JxW_values - | update_JxW_values; + out |= update_gradients | update_piola | update_covariant_transformation; + if (flags & update_hessians) - out |= update_hessians - | update_contravariant_transformation - | update_covariant_transformation; - + out |= update_hessians | update_piola | update_covariant_transformation; + return out; } diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc index 8af15830ad..9c0087f7cd 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc @@ -42,7 +42,7 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg) Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_none; + this->mapping_type = mapping_raviart_thomas; // These must be done first, since // they change the evaluation of // basis functions @@ -770,6 +770,33 @@ FE_RaviartThomasNodal<dim>::initialize_support_points (const unsigned int deg) } +template <int dim> +UpdateFlags +FE_RaviartThomasNodal<dim>::update_once (const UpdateFlags) const +{ + return update_default; +} + + +template <int dim> +UpdateFlags +FE_RaviartThomasNodal<dim>::update_each (const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + if (flags & update_values) + out |= update_values | update_piola; + + if (flags & update_gradients) + out |= update_gradients | update_piola | update_covariant_transformation; + + if (flags & update_hessians) + out |= update_hessians | update_piola | update_covariant_transformation; + + return out; +} + + template class FE_RaviartThomasNodal<deal_II_dimension>; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/fe/mapping.cc b/deal.II/deal.II/source/fe/mapping.cc index 5a2ae42186..e0d8290ede 100644 --- a/deal.II/deal.II/source/fe/mapping.cc +++ b/deal.II/deal.II/source/fe/mapping.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2005, 2006 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2005, 2006, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -22,6 +22,69 @@ Mapping<dim, spacedim>::~Mapping () {} +// This function is deprecated and has been replaced by transform above +template<int dim, int spacedim> +void +Mapping<dim,spacedim>::transform_covariant ( + const VectorSlice<const std::vector<Tensor<1,dim> > > input, + const unsigned int offset, + VectorSlice<std::vector<Tensor<1,spacedim> > > output, + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const +{ + Assert (offset == 0, ExcInternalError()); + + transform(input, output, mapping_data, mapping_covariant); +} + + + +// This function is deprecated and has been replaced by transform above +template <int dim, int spacedim> +void +Mapping<dim, spacedim>::transform_covariant ( + const VectorSlice<const std::vector<Tensor<2,dim> > > input, + const unsigned int offset, + VectorSlice<std::vector<Tensor<2,spacedim> > > output, + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const +{ + Assert (offset == 0, ExcInternalError()); + + transform(input, output, mapping_data, mapping_covariant); +} + + + +// This function is deprecated and has been replaced by transform above +template<int dim, int spacedim> +void +Mapping<dim,spacedim>::transform_contravariant ( + const VectorSlice<const std::vector<Tensor<1,dim> > > input, + const unsigned int offset, + VectorSlice<std::vector<Tensor<1,spacedim> > > output, + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const +{ + Assert (offset == 0, ExcInternalError()); + + transform(input, output, mapping_data, mapping_contravariant); +} + + + +// This function is deprecated and has been replaced by transform above +template<int dim, int spacedim> +void +Mapping<dim,spacedim>::transform_contravariant ( + const VectorSlice<const std::vector<Tensor<2,dim> > > input, + const unsigned int offset, + VectorSlice<std::vector<Tensor<2,spacedim> > > output, + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const +{ + Assert (offset == 0, ExcInternalError()); + + transform(input, output, mapping_data, mapping_contravariant); +} + + /*------------------------------ InternalDataBase ------------------------------*/ diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index c9f4dbb9c3..20fad4a21d 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -534,40 +534,38 @@ MappingCartesian<dim,spacedim>::transform ( switch (mapping_type) { case mapping_covariant: - if (true) - { - Assert (data.update_flags & update_covariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - for (unsigned int i=0; i<output.size(); ++i) - for (unsigned int d=0;d<dim;++d) - output[i][d] = input[i][d]/data.length[d]; - return; - } + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i<output.size(); ++i) + for (unsigned int d=0;d<dim;++d) + output[i][d] = input[i][d]/data.length[d]; + return; + } + case mapping_contravariant: - if (true) - { - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - for (unsigned int i=0; i<output.size(); ++i) - for (unsigned int d=0;d<dim;++d) - output[i][d] = input[i][d]*data.length[d]; - return; - } + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i<output.size(); ++i) + for (unsigned int d=0;d<dim;++d) + output[i][d] = input[i][d]*data.length[d]; + return; + } case mapping_piola: - if (true) - { - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - Assert (data.update_flags & update_volume_elements, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - for (unsigned int i=0; i<output.size(); ++i) - for (unsigned int d=0;d<dim;++d) - output[i][d] = input[i][d] * data.length[d] / data.volume_element; - return; - } + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i<output.size(); ++i) + for (unsigned int d=0;d<dim;++d) + output[i][d] = input[i][d] * data.length[d] / data.volume_element; + return; + } default: Assert(false, ExcNotImplemented()); } @@ -578,136 +576,82 @@ void MappingCartesian<dim,spacedim>::transform ( const VectorSlice<const std::vector<Tensor<2,dim> > > input, VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal, + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data, const MappingType mapping_type) const { + AssertDimension (input.size(), output.size()); + Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, + ExcInternalError()); + const InternalData &data = static_cast<const InternalData&>(mapping_data); + switch (mapping_type) { case mapping_covariant: - transform_covariant(input, 0, output, internal); - return; - case mapping_contravariant: - transform_contravariant(input, 0, output, internal); - return; - default: - Assert(false, ExcNotImplemented()); - } -} - - - -template<int dim, int spacedim> -void -MappingCartesian<dim, spacedim>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - const InternalData &data = dynamic_cast<const InternalData&> (mapping_data); - - Assert (data.update_flags & update_covariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - // simply scale by inverse Jacobian - // (which is diagonal here) - for (unsigned int i=0; i<output.size(); ++i) - for (unsigned int d=0;d<dim;++d) - output[i][d] = input[i+offset][d]/data.length[d]; -} - - - -template<int dim, int spacedim> -void -MappingCartesian<dim, spacedim>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - const InternalData &data = dynamic_cast<const InternalData&> (mapping_data); - - Assert (data.update_flags & update_covariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - Assert (output.size() + offset <= input.size(), ExcInternalError()); - - // simply scale by inverse Jacobian - // (which is diagonal here) - for (unsigned int i=0; i<output.size(); ++i) - for (unsigned int d=0; d<dim; ++d) - for (unsigned int p=0; p<dim; ++p) - output[i][d][p] = input[i+offset][d][p]/data.length[p]; -} - - - -template<int dim, int spacedim> -void -MappingCartesian<dim, spacedim>::transform_contravariant ( - const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - const InternalData &data = dynamic_cast<const InternalData&> (mapping_data); - - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - Assert (output.size() + offset <= input.size(), ExcInternalError()); - - // simply scale by inverse Jacobian - // (which is diagonal here) - for (unsigned int i=0; i<output.size(); ++i) - for (unsigned int d=0;d<dim;++d) - output[i][d] = input[i+offset][d] * data.length[d]; -} - - + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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()); + + 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; + } -template<int dim, int spacedim> -void -MappingCartesian<dim, spacedim>::transform_contravariant ( - const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - const InternalData &data = dynamic_cast<const InternalData&> (mapping_data); + case mapping_covariant_gradient: + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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; + } - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + case mapping_contravariant_gradient: + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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; + } - Assert (output.size() + offset <= input.size(), ExcInternalError()); - - // simply scale by inverse Jacobian - // (which is diagonal here) - for (unsigned int i=0; i<output.size(); ++i) - for (unsigned int d=0; d<dim; ++d) - for (unsigned int p=0; p<dim; ++p) - output[i][d][p] = input[i+offset][d][p] * data.length[p]; + case mapping_piola: + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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] / data.volume_element; + return; + } + + default: + Assert(false, ExcNotImplemented()); + } } diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index 53e84d17b0..22940c6486 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -1207,198 +1207,74 @@ void MappingQ<dim,spacedim>::transform ( const VectorSlice<const std::vector<Tensor<1,dim> > > input, VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal, + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data, const MappingType mapping_type) const { - switch (mapping_type) - { - case mapping_covariant: - transform_covariant(input, 0, output, internal); - return; - case mapping_contravariant: - transform_contravariant(input, 0, output, internal); - return; - default: - Assert(false, ExcNotImplemented()); - } -} - - -template<int dim, int spacedim> -void -MappingQ<dim,spacedim>::transform ( - const VectorSlice<const std::vector<Tensor<2,dim> > > input, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal, - const MappingType mapping_type) const -{ - switch (mapping_type) - { - case mapping_covariant: - transform_covariant(input, 0, output, internal); - return; - case mapping_contravariant: - transform_contravariant(input, 0, output, internal); - return; - default: - Assert(false, ExcNotImplemented()); - } -} - - - -template<int dim, int spacedim> -void -MappingQ<dim,spacedim>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<1,spacedim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - const typename MappingQ1<dim,spacedim>::InternalData *q1_data = - dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data); - Assert(q1_data!=0, ExcInternalError()); - - typename std::vector<Tensor<2,dim> >::const_iterator tensor; - - if (q1_data->is_mapping_q1_data) - tensor = q1_data->covariant.begin(); - else - { - const InternalData *data = dynamic_cast<const InternalData *> (q1_data); - Assert(data!=0, ExcInternalError()); - - if (data->use_mapping_q1_on_current_cell) - tensor = data->mapping_q1_data.covariant.begin(); - else - tensor = data->covariant.begin(); - } - - Tensor<1, spacedim> auxiliary; - - for (unsigned int i=0; i<output.size(); ++i) - { - for (unsigned int d=0;d<dim;++d) - auxiliary[d] = input[i][d]; - contract (output[i], auxiliary, *(tensor++)); - } -} - - - -template<int dim, int spacedim> -void -MappingQ<dim,spacedim>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<2,spacedim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - const typename MappingQ1<dim,spacedim>::InternalData *q1_data = - dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data); - Assert(q1_data!=0, ExcInternalError()); - - typename std::vector<Tensor<2,dim> >::const_iterator tensor; - - if (q1_data->is_mapping_q1_data) - tensor = q1_data->covariant.begin(); - else - { - const InternalData *data = dynamic_cast<const InternalData *> (q1_data); - Assert(data!=0, ExcInternalError()); - - if (data->use_mapping_q1_on_current_cell) - tensor = data->mapping_q1_data.covariant.begin(); - else - tensor = data->covariant.begin(); - } - - for (unsigned int i=0; i<output.size(); ++i) - contract (output[i], input[i+offset], *(tensor++)); -} - - - -template<int dim, int spacedim> -void -MappingQ<dim,spacedim>::transform_contravariant ( - const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); AssertDimension (input.size(), output.size()); - + // The data object may be jsut a + // MappingQ1::InternalData, so we + // have to test for this first. const typename MappingQ1<dim,spacedim>::InternalData *q1_data = dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data); Assert(q1_data!=0, ExcInternalError()); - - typename std::vector<Tensor<2,dim> >::const_iterator tensor; - if (q1_data->is_mapping_q1_data) - tensor = q1_data->contravariant.begin(); - else - { - const InternalData *data = dynamic_cast<const InternalData *> (q1_data); - Assert(data!=0, ExcInternalError()); - - if (data->use_mapping_q1_on_current_cell) - tensor = data->mapping_q1_data.contravariant.begin(); - else - tensor = data->contravariant.begin(); - } - - Tensor<1, spacedim> auxiliary; - - for (unsigned int i=0; i<output.size(); ++i) + // If it is a genuine + // MappingQ::InternalData, we have + // to test further + if (!q1_data->is_mapping_q1_data) { - for (unsigned int d=0;d<dim;++d) - auxiliary[d] = input[i][d]; - contract (output[i], *(tensor++), auxiliary); + Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, + ExcInternalError()); + const InternalData &data = static_cast<const InternalData&>(mapping_data); + // If we only use the + // Q1-portion, we have to + // extract that data object + if (data.use_mapping_q1_on_current_cell) + q1_data = &data.mapping_q1_data; } + // Now, q1_data should have the + // right tensors in it and we call + // the base classes transform + // function + MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type); } - template<int dim, int spacedim> void -MappingQ<dim,spacedim>::transform_contravariant ( +MappingQ<dim,spacedim>::transform ( const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data, + const MappingType mapping_type) const { - Assert (offset == 0, ExcInternalError()); AssertDimension (input.size(), output.size()); - + // The data object may be jsut a + // MappingQ1::InternalData, so we + // have to test for this first. const typename MappingQ1<dim,spacedim>::InternalData *q1_data = dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data); Assert(q1_data!=0, ExcInternalError()); - - typename std::vector<Tensor<2,dim> >::const_iterator tensor; - if (q1_data->is_mapping_q1_data) - tensor = q1_data->contravariant.begin(); - else + // If it is a genuine + // MappingQ::InternalData, we have + // to test further + if (!q1_data->is_mapping_q1_data) { - const InternalData *data = dynamic_cast<const InternalData *> (q1_data); - Assert(data!=0, ExcInternalError()); - - if (data->use_mapping_q1_on_current_cell) - tensor = data->mapping_q1_data.contravariant.begin(); - else - tensor = data->contravariant.begin(); + Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, + ExcInternalError()); + const InternalData &data = static_cast<const InternalData&>(mapping_data); + // If we only use the + // Q1-portion, we have to + // extract that data object + if (data.use_mapping_q1_on_current_cell) + q1_data = &data.mapping_q1_data; } - - for (unsigned int i=0; i<output.size(); ++i) - contract (output[i], *(tensor++), input[i+offset]); + // Now, q1_data should have the + // right tensors in it and we call + // the base classes transform + // function + MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type); } diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index d229bda802..740274e8c6 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -729,15 +729,16 @@ MappingQ1<dim,spacedim>::compute_mapping_support_points( template<int dim, int spacedim> void -MappingQ1<dim,spacedim>::fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell, - const Quadrature<dim> &q, - typename Mapping<dim,spacedim>::InternalDataBase &mapping_data, - std::vector<Point<spacedim> > &quadrature_points, - std::vector<double> &JxW_values, - std::vector<Tensor<2,spacedim> > &jacobians, - std::vector<Tensor<3,spacedim> > &jacobian_grads, - std::vector<Tensor<2,spacedim> > &inverse_jacobians, - std::vector<Point<spacedim> > &cell_normal_vectors) const +MappingQ1<dim,spacedim>::fill_fe_values ( + const typename Triangulation<dim,spacedim>::cell_iterator &cell, + const Quadrature<dim> &q, + typename Mapping<dim,spacedim>::InternalDataBase &mapping_data, + std::vector<Point<spacedim> > &quadrature_points, + std::vector<double> &JxW_values, + std::vector<Tensor<2,spacedim> > &jacobians, + std::vector<Tensor<3,spacedim> > &jacobian_grads, + std::vector<Tensor<2,spacedim> > &inverse_jacobians, + std::vector<Point<spacedim> > &cell_normal_vectors) const { // ensure that the following cast is really correct: Assert (dynamic_cast<InternalData *>(&mapping_data) != 0, @@ -1184,49 +1185,49 @@ MappingQ1<dim,spacedim>::transform ( switch (mapping_type) { case mapping_covariant: - if (true) - { - Assert (data.update_flags & update_covariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - for (unsigned int i=0; i<output.size(); ++i) - { - for (unsigned int d=0;d<dim;++d) - auxiliary[d] = input[i][d]; - contract (output[i], auxiliary, data.covariant[i]); - } - return; - } + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i<output.size(); ++i) + { + for (unsigned int d=0;d<dim;++d) + auxiliary[d] = input[i][d]; + contract (output[i], auxiliary, data.covariant[i]); + } + return; + } + case mapping_contravariant: - if (true) - { - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - for (unsigned int i=0; i<output.size(); ++i) - { - for (unsigned int d=0;d<dim;++d) - auxiliary[d] = input[i][d]; - contract (output[i], data.contravariant[i], auxiliary); - } - return; - } + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i<output.size(); ++i) + { + for (unsigned int d=0;d<dim;++d) + auxiliary[d] = input[i][d]; + contract (output[i], data.contravariant[i], auxiliary); + } + return; + } + case mapping_piola: - if (true) - { - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - Assert (data.update_flags & update_volume_elements, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - for (unsigned int i=0; i<output.size(); ++i) - { - for (unsigned int d=0;d<dim;++d) - auxiliary[d] = input[i][d] / data.volume_elements[i]; - contract (output[i], data.contravariant[i], auxiliary); - } - return; - } + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i<output.size(); ++i) + { + for (unsigned int d=0;d<dim;++d) + auxiliary[d] = input[i][d] / data.volume_elements[i]; + contract (output[i], data.contravariant[i], auxiliary); + } + return; + } + default: Assert(false, ExcNotImplemented()); } @@ -1238,206 +1239,114 @@ void MappingQ1<dim,spacedim>::transform ( const VectorSlice<const std::vector<Tensor<2, dim> > > input, VectorSlice<std::vector<Tensor<2, spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &internal, + const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data, const MappingType mapping_type) const { AssertDimension (input.size(), output.size()); - - switch (mapping_type) - { - case mapping_covariant: - transform_covariant(input, 0, output, internal); - return; - case mapping_contravariant: - transform_contravariant(input, 0, output, internal); - return; - default: - Assert(false, ExcNotImplemented()); - } -} - - -// This function is deprecated and has been replaced by transform above -template<int dim, int spacedim> -void -MappingQ1<dim,spacedim>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - // ensure that the following cast is really correct: - Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast<const InternalData&>(mapping_data); - - Assert (data.update_flags & update_covariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - Tensor<1, spacedim> auxiliary; - - for (unsigned int i=0; i<output.size(); ++i) - { - for (unsigned int d=0;d<dim;++d) - auxiliary[d] = input[i][d]; - contract (output[i], auxiliary, data.covariant[i]); - } -} - - - -// This function is deprecated and has been replaced by transform above -template <int dim, int spacedim> -void -MappingQ1<dim, spacedim>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const -{ - Assert (offset == 0, ExcInternalError()); - AssertDimension (input.size(), output.size()); - - // ensure that the following cast - // is really correct: - Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast<const InternalData&>(mapping_data); - - Assert (data.update_flags & update_covariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - for (unsigned int i=0; i<output.size(); ++i) - contract (output[i], input[i+offset], data.covariant[i]); -} - -#if deal_II_dimension == 1 - - -// This function is deprecated and has been replaced by transform above -template <> -void -MappingQ1<1, 2>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<2,1> > > , - const unsigned int, - VectorSlice<std::vector<Tensor<2,2> > > , - const Mapping<1,2>::InternalDataBase &) const -{ - Assert(false, ExcNotImplemented()); -} - -#endif - -#if deal_II_dimension == 2 - - -// This function is deprecated and has been replaced by transform above -template <> -void -MappingQ1<2, 3>::transform_covariant ( - const VectorSlice<const std::vector<Tensor<2,2> > > , - const unsigned int, - VectorSlice<std::vector<Tensor<2,3> > > , - const Mapping<2,3>::InternalDataBase &) const -{ - Assert(false, ExcNotImplemented()); -} - -#endif - - -// This function is deprecated and has been replaced by transform above -template<int dim, int spacedim> -void -MappingQ1<dim,spacedim>::transform_contravariant ( - const VectorSlice<const std::vector<Tensor<1,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<1,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const -{ - // ensure that the following cast is really correct: Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, ExcInternalError()); const InternalData &data = static_cast<const InternalData&>(mapping_data); - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - Assert (output.size() + offset <= input.size(), ExcInternalError()); - - Tensor<1, spacedim> auxiliary; + Tensor<2, spacedim> aux1; + Tensor<2, spacedim> aux2; - for (unsigned int i=0; i<output.size(); ++i) + switch (mapping_type) { - for (unsigned int d=0;d<dim;++d) - auxiliary[d] = input[i][d]; - contract (output[i], data.contravariant[i], auxiliary); + case mapping_covariant: + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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) + aux1[d1][d2] = input[i][d1][d2]; + contract (output[i], aux1, data.covariant[i]); + } + return; + } + + case mapping_contravariant: + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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) + aux1[d1][d2] = input[i][d1][d2]; + contract (output[i], data.contravariant[i], aux1); + } + return; + } + + case mapping_covariant_gradient: + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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) + aux1[d1][d2] = input[i][d1][d2]; + + contract(aux2, aux1, data.covariant[i]); + contract(output[i], data.covariant[i], aux2); + } + + return; + } + + case mapping_contravariant_gradient: + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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) + aux1[d1][d2] = input[i][d1][d2]; + + contract(aux2, aux1, data.covariant[i]); + contract(output[i], data.contravariant[i], aux2); + } + + return; + } + + case mapping_piola: + { + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase<dim>::ExcAccessToUninitializedField()); + + 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) + aux1[d1][d2] = input[i][d1][d2] / data.volume_elements[i]; + + contract (aux2, aux1, data.covariant[i]); + contract (output[i], data.contravariant[i], aux2); + } + return; + } + default: + Assert(false, ExcNotImplemented()); } } -// This function is deprecated and has been replaced by transform above -template<int dim, int spacedim> -void -MappingQ1<dim,spacedim>::transform_contravariant ( - const VectorSlice<const std::vector<Tensor<2,dim> > > input, - const unsigned int offset, - VectorSlice<std::vector<Tensor<2,spacedim> > > output, - const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data) const -{ - // ensure that the following cast is really correct: - Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast<const InternalData&>(mapping_data); - - Assert (data.update_flags & update_contravariant_transformation, - typename FEValuesBase<dim>::ExcAccessToUninitializedField()); - - Assert (output.size() + offset <= input.size(), ExcInternalError()); - - for (unsigned int i=0; i<output.size(); ++i) - contract (output[i], data.contravariant[i], input[i+offset]); -} - - -#if deal_II_dimension == 1 - -// This function is deprecated and has been replaced by transform above -template <> -void -MappingQ1<1, 2>::transform_contravariant ( - const VectorSlice<const std::vector<Tensor<2,1> > > , - const unsigned int, - VectorSlice<std::vector<Tensor<2,2> > > , - const Mapping<1,2>::InternalDataBase &) const -{ - Assert(false, ExcNotImplemented()); -} - -#endif - -#if deal_II_dimension == 2 - -// This function is deprecated and has been replaced by transform above -template <> -void -MappingQ1<2, 3>::transform_contravariant ( - const VectorSlice<const std::vector<Tensor<2,2> > > , - const unsigned int, - VectorSlice<std::vector<Tensor<2,3> > > , - const Mapping<2,3>::InternalDataBase &) const -{ - Assert(false, ExcNotImplemented()); -} - -#endif - - - template<int dim, int spacedim> Point<spacedim> MappingQ1<dim,spacedim>::transform_unit_to_real_cell ( -- 2.39.5