From: wolf Date: Tue, 3 Sep 2002 15:30:31 +0000 (+0000) Subject: Add and clarify some comments. Add a second set of transform_{co,contra}variant funct... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7088512d3d0ba3cb57599f09ee835d7d29ed9c1e;p=dealii-svn.git Add and clarify some comments. Add a second set of transform_{co,contra}variant functions that transform matrices instead of vectors. git-svn-id: https://svn.dealii.org/trunk@6349 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index e0213bc5f3..52618daa7f 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -18,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -195,36 +194,97 @@ class Mapping : public Subscriptor }; /** - * Tranform a field of covariant - * vectors. The vector spanned by - * @p{begin} and @p{end} must + * Transform a field of covariant + * vectors. The covariant + * transformation multiplies a + * vector from the right with the + * inverse of the Jacobian matrix + * of the transformation from + * unit to real space + * cell. Alternatively, this is + * equivalent to a multiplication + * from the left with the + * transpose if the inverse + * matrix. + * + * The range of vectors spanned + * by @p{begin} and @p{end} must * have as many elements as there * are quadrature points (not * tested inside the function). + * Also note the different + * convention for parameters + * compared to the standard C++ + * @p{transform} function: here, + * first destination, then source + * are specified, and the number + * of elements is determined by a + * range of destination + * vectors. This convention is + * due to the way the function is + * usually called. * * The vector @p{src} must * contain at least as many - * elements. + * elements as there are + * quadrature points. */ virtual void transform_covariant (Tensor<1,dim> *begin, Tensor<1,dim> *end, const Tensor<1,dim> *src, + const InternalDataBase &internal) const = 0; + + /** + * Transform a set of matrices + * covariantly, i.e. multiply + * each function in the input + * range by the Jacobian matrices + * at the different quadrature + * points from the left. + * + * The same applies as to the + * function above regarding input + * and output ranges. + */ + virtual + void + transform_covariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, const InternalDataBase &internal) const = 0; /** - * Tranform a field of + * Transform a field of * contravariant vectors. The - * vector spanned by @p{begin} - * and @p{end} must have as many - * elements as there are - * quadrature points (not tested - * inside the function). + * contravariant transformation + * multiplies a vector from the + * left with the Jacobian matrix + * of the transformation from + * unit to real space cell. + * + * The range of vectors spanned + * by @p{begin} and @p{end} must + * have as many elements as there + * are quadrature points (not + * tested inside the function). + * Also note the different + * convention for parameters + * compared to the standard C++ + * @p{transform} function: here, + * first destination, then source + * are specified, and the number + * of elements is determined by a + * range of destination + * vectors. This convention is + * due to the way the function is + * usually called. * * The vector @p{src} must * contain at least as many - * elements. + * elements as there are + * quadrature points. */ virtual void @@ -233,6 +293,32 @@ class Mapping : public Subscriptor const Tensor<1,dim> *src, const InternalDataBase &internal) const = 0; + /** + * Transform a set of matrices + * contravariantly, i.e. multiply + * each function in the input + * range by the inverse Jacobian + * matrices at the different + * quadrature points from the + * right. Note that here it is no + * more equivalent to + * multiplication with the + * transpose of the inverse + * matrices from the left, unless + * the matrices to be multiplied + * with are symmetric. + * + * The same applies as to the + * function above regarding input + * and output ranges. + */ + virtual + void + transform_contravariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const InternalDataBase &internal) const = 0; + /** * Indicate fields to be updated * in the constructor of @@ -352,7 +438,7 @@ class Mapping : public Subscriptor fill_fe_values (const typename DoFHandler::cell_iterator &cell, const Quadrature &quadrature, InternalDataBase &internal, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values) const = 0; /** diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index 0027696eb5..134943693b 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -108,6 +108,16 @@ class MappingCartesian : public Mapping transform_covariant (Tensor<1,dim> *begin, Tensor<1,dim> *end, const Tensor<1,dim> *src, + const typename Mapping::InternalDataBase &internal) const; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + transform_covariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, const typename Mapping::InternalDataBase &internal) const; /** @@ -118,6 +128,16 @@ class MappingCartesian : public Mapping transform_contravariant (Tensor<1,dim> *begin, Tensor<1,dim> *end, const Tensor<1,dim> *src, + const typename Mapping::InternalDataBase &internal) const; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + transform_contravariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, const typename Mapping::InternalDataBase &internal) const; /** @@ -205,12 +225,12 @@ class MappingCartesian : public Mapping * * Filled once. */ - vector2d > unit_tangentials; + Table<2,Tensor<1,dim> > unit_tangentials; /** * Auxiliary vectors for internal use. */ - vector2d > aux; + Table<2,Tensor<1,dim> > aux; }; /** diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index a935afb806..dbc26e4a1e 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -15,7 +15,7 @@ #include -#include +#include #include template class TensorProductPolynomials; @@ -79,8 +79,8 @@ class MappingQ : public MappingQ1 const Point &p) const; /** - * Implementation of the interface in - * @ref{Mapping}. + * Implementation of the + * interface in @ref{Mapping}. */ virtual void transform_covariant (Tensor<1,dim> *begin, @@ -89,8 +89,18 @@ class MappingQ : public MappingQ1 const typename Mapping::InternalDataBase &internal) const; /** - * Implementation of the interface in - * @ref{Mapping}. + * Implementation of the + * interface in @ref{Mapping}. + */ + virtual void + transform_covariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &internal) const; + + /** + * Implementation of the + * interface in @ref{Mapping}. */ virtual void transform_contravariant (Tensor<1,dim> *begin, @@ -98,6 +108,16 @@ class MappingQ : public MappingQ1 const Tensor<1,dim> *src, const typename Mapping::InternalDataBase &internal) const; + /** + * Implementation of the + * interface in @ref{Mapping}. + */ + virtual void + transform_contravariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &internal) const; + /** * Return the degree of the * mapping, i.e. the value which @@ -304,7 +324,7 @@ class MappingQ : public MappingQ1 * this vector is computed. */ void - set_laplace_on_quad_vector(vector2d &loqvs) const; + set_laplace_on_quad_vector(Table<2,double> &loqvs) const; /** * This function is needed by the @@ -317,7 +337,7 @@ class MappingQ : public MappingQ1 * @p{degree>2} this vector is * computed. */ - void set_laplace_on_hex_vector(vector2d &lohvs) const; + void set_laplace_on_hex_vector(Table<2,double> &lohvs) const; /** * Computes the @p{laplace_on_quad(hex)_vector}. @@ -327,7 +347,7 @@ class MappingQ : public MappingQ1 * functions if the data is not * yet hardcoded. */ - void compute_laplace_vector(vector2d &lvs) const; + void compute_laplace_vector(Table<2,double> &lvs) const; /** * Takes a @@ -344,7 +364,7 @@ class MappingQ : public MappingQ1 * @p{n_inner} computed inner * points are appended. */ - void apply_laplace_vector(const vector2d &lvs, + void apply_laplace_vector(const Table<2,double> &lvs, std::vector > &a) const; /** @@ -426,7 +446,7 @@ class MappingQ : public MappingQ1 * unit_support_points on the * boundary of the quad */ - vector2d laplace_on_quad_vector; + Table<2,double> laplace_on_quad_vector; /** * Needed by the @@ -434,7 +454,7 @@ class MappingQ : public MappingQ1 * (for @p{dim==3}). Filled by the * constructor. */ - vector2d laplace_on_hex_vector; + Table<2,double> laplace_on_hex_vector; /** * Exception. @@ -514,11 +534,11 @@ template<> void MappingQ<1>::compute_shapes_virtual ( const std::vector > &unit_points, MappingQ1<1>::InternalData &data) const; template <> void MappingQ<1>::set_laplace_on_quad_vector( - vector2d &) const; + Table<2,double> &) const; template <> void MappingQ<3>::set_laplace_on_hex_vector( - vector2d &lohvs) const; + Table<2,double> &lohvs) const; template <> void MappingQ<1>::compute_laplace_vector( - vector2d &) const; + Table<2,double> &) const; template <> void MappingQ<1>::add_line_support_points ( const Triangulation<1>::cell_iterator &, std::vector > &) const; diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index 24237bb452..6aa9655843 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -15,6 +15,7 @@ #include +#include #include #include @@ -82,9 +83,29 @@ class MappingQ1 : public Mapping * @ref{Mapping}. */ virtual void + transform_covariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &internal) const; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void transform_contravariant (Tensor<1,dim> *begin, Tensor<1,dim> *end, const Tensor<1,dim> *src, + const typename Mapping::InternalDataBase &internal) const; + + /** + * Implementation of the interface in + * @ref{Mapping}. + */ + virtual void + transform_contravariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, const typename Mapping::InternalDataBase &internal) const; /** @@ -237,7 +258,14 @@ class MappingQ1 : public Mapping /** * Tensors of covariant - * transformation. + * transformation at each of + * the quadrature points. The + * matrix stored is the + * inverse of the Jacobian + * matrix, which itself is + * stored in the + * @p{contravariant} field of + * this structure. * * Computed on each cell. */ @@ -245,7 +273,12 @@ class MappingQ1 : public Mapping /** * Tensors of covariant - * transformation. + * transformation at each of + * the quadrature points. The + * contravariant matrix is + * the Jacobian of the + * transformation, + * i.e. $J_ij=dx_i/d\hat x_j$. * * Computed on each cell. */ @@ -259,12 +292,12 @@ class MappingQ1 : public Mapping * * Filled once. */ - vector2d > unit_tangentials; + Table<2,Tensor<1,dim> > unit_tangentials; /** * Auxuliary vectors for internal use. */ - vector2d > aux; + Table<2,Tensor<1,dim> > aux; /** * Stores the support points of diff --git a/deal.II/deal.II/include/fe/mapping_q1_eulerian.h b/deal.II/deal.II/include/fe/mapping_q1_eulerian.h index 003cc10615..c2d0260095 100644 --- a/deal.II/deal.II/include/fe/mapping_q1_eulerian.h +++ b/deal.II/deal.II/include/fe/mapping_q1_eulerian.h @@ -95,8 +95,8 @@ class MappingQ1Eulerian : public MappingQ1 * can be initialized by * @p{DoFObjectAccessor::set_dof_values()}. */ - MappingQ1Eulerian ( const Vector &euler_transform_vectors, - const DoFHandler &shiftmap_dof_handler); + MappingQ1Eulerian (const Vector &euler_transform_vectors, + const DoFHandler &shiftmap_dof_handler); /** diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index cc90555189..7fa2d74579 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -460,6 +460,33 @@ MappingCartesian::transform_covariant (Tensor<1,dim> *begin, } + +template +void +MappingCartesian::transform_covariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &mapping_data) const +{ + const InternalData &data = dynamic_cast (mapping_data); + + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField()); + + // simply scale by inverse Jacobian + // (which is diagonal here) + while (begin!=end) + { + for (unsigned int d=0; d void MappingCartesian::transform_contravariant (Tensor<1,dim> *begin, @@ -488,6 +515,37 @@ MappingCartesian::transform_contravariant (Tensor<1,dim> *begin, } + +template +void +MappingCartesian::transform_contravariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &mapping_data) const +{ + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + const InternalData &data = dynamic_cast (mapping_data); + + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField()); + + // simply scale by Jacobian + // (which is diagonal here) + while (begin!=end) + { + for (unsigned int d=0; d Point MappingCartesian::transform_unit_to_real_cell ( const typename Triangulation::cell_iterator cell, diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index d55bee4885..b0d5ac2ab7 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -61,8 +61,6 @@ MappingQ::InternalData::memory_consumption () const // cells are scaled linearly template<> MappingQ<1>::MappingQ (const unsigned int): - laplace_on_quad_vector(0), - laplace_on_hex_vector(0), degree(1), n_inner(0), n_outer(0), @@ -100,9 +98,8 @@ static number power(const number x, const unsigned int y) template -MappingQ::MappingQ (const unsigned int p): - laplace_on_quad_vector(0), - laplace_on_hex_vector(0), +MappingQ::MappingQ (const unsigned int p) + : degree(p), n_inner(power(degree-1, dim)), n_outer((dim==2) ? 4+4*(degree-1) @@ -409,7 +406,7 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter template <> void -MappingQ<1>::set_laplace_on_quad_vector(vector2d &) const +MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const { Assert(false, ExcInternalError()); } @@ -418,7 +415,7 @@ MappingQ<1>::set_laplace_on_quad_vector(vector2d &) const template void -MappingQ::set_laplace_on_quad_vector(vector2d &loqvs) const +MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const { Assert(degree>1, ExcInternalError()); const unsigned int n_inner_2d=(degree-1)*(degree-1); @@ -490,7 +487,7 @@ MappingQ::set_laplace_on_quad_vector(vector2d &loqvs) const template <> void -MappingQ<3>::set_laplace_on_hex_vector(vector2d &lohvs) const +MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const { Assert(degree>1, ExcInternalError()); @@ -536,7 +533,7 @@ MappingQ<3>::set_laplace_on_hex_vector(vector2d &lohvs) const template void -MappingQ::set_laplace_on_hex_vector(vector2d &) const +MappingQ::set_laplace_on_hex_vector(Table<2,double> &) const { Assert(false, ExcInternalError()); } @@ -548,7 +545,7 @@ MappingQ::set_laplace_on_hex_vector(vector2d &) const template <> void -MappingQ<1>::compute_laplace_vector(vector2d &) const +MappingQ<1>::compute_laplace_vector(Table<2,double> &) const { Assert(false, ExcInternalError()); } @@ -558,7 +555,7 @@ MappingQ<1>::compute_laplace_vector(vector2d &) const template void -MappingQ::compute_laplace_vector(vector2d &lvs) const +MappingQ::compute_laplace_vector(Table<2,double> &lvs) const { Assert(lvs.n_rows()==0, ExcInternalError()); Assert(dim==2 || dim==3, ExcNotImplemented()); @@ -617,7 +614,7 @@ MappingQ::compute_laplace_vector(vector2d &lvs) const template void -MappingQ::apply_laplace_vector(const vector2d &lvs, +MappingQ::apply_laplace_vector(const Table<2,double> &lvs, std::vector > &a) const { Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree)); @@ -1086,6 +1083,40 @@ MappingQ::transform_covariant (Tensor<1,dim> *begin, } + +template +void +MappingQ::transform_covariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &mapping_data) const +{ + const typename MappingQ1::InternalData *q1_data = + dynamic_cast::InternalData *> (&mapping_data); + Assert(q1_data!=0, ExcInternalError()); + + typename std::vector >::const_iterator tensor; + + if (q1_data->is_mapping_q1_data) + tensor = q1_data->covariant.begin(); + else + { + const InternalData *data = dynamic_cast (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(); + } + while (begin!=end) + { + contract (*(begin++), *(src++), *(tensor++)); + } +} + + + template void MappingQ::transform_contravariant (Tensor<1,dim> *begin, @@ -1118,6 +1149,40 @@ MappingQ::transform_contravariant (Tensor<1,dim> *begin, } + +template +void +MappingQ::transform_contravariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &mapping_data) const +{ + const typename MappingQ1::InternalData *q1_data = + dynamic_cast::InternalData *> (&mapping_data); + Assert(q1_data!=0, ExcInternalError()); + + typename std::vector >::const_iterator tensor; + + if (q1_data->is_mapping_q1_data) + tensor = q1_data->contravariant.begin(); + else + { + const InternalData *data = dynamic_cast (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(); + } + while (begin!=end) + { + contract (*(begin++), *(tensor++), *(src++)); + } +} + + + template Point MappingQ::transform_unit_to_real_cell ( const typename Triangulation::cell_iterator cell, diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index bc31b61396..8a31a3129b 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -467,47 +467,45 @@ MappingQ1::compute_fill (const typename DoFHandler::cell_iterator &cel // ExcDimensionMismatch(covariant_grads.size(), npts)); } - std::vector > &a=data.mapping_support_points; - - // store all Lagrangian - // support points in a - if (a.size()==0 - || (&cell->get_triangulation() != - &data.cell_of_current_support_points->get_triangulation()) - || (cell!=data.cell_of_current_support_points)) + // if necessary, recompute the + // support points of the + // transformation of this cell + if ((data.mapping_support_points.size() == 0) + || + (&cell->get_triangulation() != + &data.cell_of_current_support_points->get_triangulation()) + || + (cell != data.cell_of_current_support_points)) { - compute_mapping_support_points(cell, a); - data.cell_of_current_support_points=cell; + compute_mapping_support_points(cell, data.mapping_support_points); + data.cell_of_current_support_points = cell; } + + // first compute quadrature points + if (update_flags & update_q_points) + for (unsigned int point=0; point::transform_covariant (Tensor<1,dim> *begin, // Assert (dst.size() == data.contravariant.size(), // ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size())); - typename std::vector >::const_iterator tensor = data.covariant.begin(); + typename std::vector >::const_iterator + tensor = data.covariant.begin(); + + while (begin!=end) + { + contract (*(begin++), *(src++), *(tensor++)); + } +} + + + +template +void +MappingQ1::transform_covariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &mapping_data) const +{ + const InternalData *data_ptr = dynamic_cast (&mapping_data); + Assert(data_ptr!=0, ExcInternalError()); + const InternalData &data=*data_ptr; + + Assert (data.update_flags & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField()); + +//TODO: [GK] Can we do a similar assertion? +// Assert (dst.size() == data.contravariant.size(), +// ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size())); + + typename std::vector >::const_iterator + tensor = data.covariant.begin(); while (begin!=end) { @@ -807,6 +835,7 @@ MappingQ1::transform_covariant (Tensor<1,dim> *begin, } + template void MappingQ1::transform_contravariant (Tensor<1,dim> *begin, @@ -821,7 +850,8 @@ MappingQ1::transform_contravariant (Tensor<1,dim> *begin, Assert (data.update_flags & update_contravariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - typename std::vector >::const_iterator tensor = data.contravariant.begin(); + typename std::vector >::const_iterator + tensor = data.contravariant.begin(); while (begin!=end) { @@ -829,6 +859,32 @@ MappingQ1::transform_contravariant (Tensor<1,dim> *begin, } } + +template +void +MappingQ1::transform_contravariant (Tensor<2,dim> *begin, + Tensor<2,dim> *end, + const Tensor<2,dim> *src, + const typename Mapping::InternalDataBase &mapping_data) const +{ + const InternalData* data_ptr = dynamic_cast (&mapping_data); + Assert(data_ptr!=0, ExcInternalError()); + const InternalData &data=*data_ptr; + + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField()); + + typename std::vector >::const_iterator + tensor = data.contravariant.begin(); + + while (begin!=end) + { + contract (*(begin++), *(tensor++), *(src++)); + } +} + + + template Point MappingQ1::transform_unit_to_real_cell ( const typename Triangulation::cell_iterator cell,