From ec2319d145e33a813081b75126eb44be611a4a81 Mon Sep 17 00:00:00 2001 From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Tue, 6 Jan 2009 21:20:12 +0000 Subject: [PATCH] use new Piola transform for vector valued functions git-svn-id: https://svn.dealii.org/trunk@18102 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_poly_tensor.h | 53 +------ deal.II/deal.II/include/fe/fe_update_flags.h | 27 ++-- deal.II/deal.II/include/fe/mapping.h | 21 ++- .../deal.II/include/fe/mapping_cartesian.h | 5 + deal.II/deal.II/source/fe/fe_abf.cc | 8 +- deal.II/deal.II/source/fe/fe_poly_tensor.cc | 131 ++++++++++-------- .../deal.II/source/fe/fe_raviart_thomas.cc | 8 +- .../source/fe/fe_raviart_thomas_nodal.cc | 4 +- .../deal.II/source/fe/mapping_cartesian.cc | 43 ++++-- deal.II/deal.II/source/fe/mapping_q1.cc | 23 +++ 10 files changed, 186 insertions(+), 137 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_poly_tensor.h b/deal.II/deal.II/include/fe/fe_poly_tensor.h index 6e44ccbd1f..e4dfabb33a 100644 --- a/deal.II/deal.II/include/fe/fe_poly_tensor.h +++ b/deal.II/deal.II/include/fe/fe_poly_tensor.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006 by the deal.II authors +// Copyright (C) 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 @@ -106,7 +106,7 @@ DEAL_II_NAMESPACE_OPEN * in #mapping_type. Therefore, each constructor should contain a line * like: * @verbatim - * this->mapping_type = this->independent_on_cartesian; + * this->mapping_type = this->mapping_none; * @endverbatim * * @see PolynomialsBDM, PolynomialsRaviartThomas @@ -219,57 +219,10 @@ class FE_PolyTensor : public FiniteElement<dim,spacedim> virtual UpdateFlags update_each (const UpdateFlags flags) const; protected: - /** - * Different options for - * transforming the basis - * functions from the reference - * cell to the actual mesh cell. - * - * Most vector valued elements - * either transform shape - * functions to keep node values - * on edges meaningful. Still, in - * special cases, it may be - * possible to avoid the mapping. - */ - enum MappingType { - /** - * No mapping has been - * selected, throw an error - * if needed. - */ - no_mapping, - /** - * Shape functions do not - * depend on actual mesh - * cell - */ - independent, - /** - * Shape functions do not - * depend on actual mesh - * cell. The mapping class - * must be - * MappingCartesian. - */ - independent_on_cartesian, - /** - * Shape functions are - * transformed covariant. - */ - covariant, - /** - * Shape functions are - * transformed - * contravariant. - */ - contravariant - }; - /** * The mapping type to be used to * map shape functions from the - * reference cell to the msh + * reference cell to the mesh * cell. */ MappingType mapping_type; diff --git a/deal.II/deal.II/include/fe/fe_update_flags.h b/deal.II/deal.II/include/fe/fe_update_flags.h index d854cf7299..b81bd855ff 100644 --- a/deal.II/deal.II/include/fe/fe_update_flags.h +++ b/deal.II/deal.II/include/fe/fe_update_flags.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -221,7 +221,12 @@ enum UpdateFlags * the Mapping. */ update_transformation_gradients = 0x2000, - + //! Determinant of the Jacobian + /** + * Compute the volume element + * in each quadrature point. + */ + update_volume_elements = 0x4000, /** * Compute the JxW values * on faces for the cell mapping @@ -231,25 +236,25 @@ enum UpdateFlags * is used in conjunction with * H_div subspaces like RT and ABF. */ - update_cell_JxW_values = 0x4000, + update_cell_JxW_values = 0x8000, /** * Update the location of the * mapped generalized support * points of the element. */ - update_support_points = 0x8000, + update_support_points = 0x10000, /** * Update the Jacobian of the * mapping in generalized * support points. */ - update_support_jacobians = 0x10000, + update_support_jacobians = 0x20000, /** * Update the inverse Jacobian * of the mapping in * generalized support points. */ - update_support_inverse_jacobians = 0x20000, + update_support_inverse_jacobians = 0x40000, /** * @deprecated Update * quadrature points @@ -266,8 +271,14 @@ enum UpdateFlags * vectors. Use * update_face_normal_vectors */ - update_normal_vectors = update_face_normal_vectors - + update_normal_vectors = update_face_normal_vectors, + //! Values needed for Piola transform + /** + * Combination of the flags + * needed for Piola transform + * of Hdiv elements. + */ + update_piola = update_volume_elements | update_contravariant_transformation }; diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index 95779156ee..51023ab054 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -32,8 +32,6 @@ template <int dim, int spacedim> class FEValues; template <int dim, int spacedim> class FEFaceValues; template <int dim, int spacedim> class FESubfaceValues; -//TODO: Offset in transform functions should be replaced by initializing VectorSlice correctly - /** * The transformation type used * for the Mapping::transform() functions. @@ -55,18 +53,21 @@ template <int dim, int spacedim> class FESubfaceValues; * for vectors. If such a * MappingType is applied to a * rank 2 tensor, it is implied - * that the resulting Tensor - * corresponds to the derivative - * of the vector. + * that the mapping is applied to + * each column. */ enum MappingType { +/// No mapping + mapping_none = 0x0000, /// Covariant mapping mapping_covariant = 0x0001, /// Contravariant mapping mapping_contravariant = 0x0002, /// The Piola transform usually used for Hdiv elements - mapping_piola = 0x0003, + mapping_piola = 0x0100, +/// The Piola transform for the derivative of an Hdiv element + mapping_piola_gradient = 0x0101, /// The mapping used for Raviart-Thomas elements mapping_raviart_thomas = mapping_piola, /// The mapping used for BDM elements @@ -261,6 +262,14 @@ class Mapping : public Subscriptor * object. */ virtual unsigned int memory_consumption () const; + + /** + * The determinant of the + * Jacobian in each + * quadrature point. Filled + * if #update_volume_elements. + */ + std::vector<double> volume_elements; /** * The positions of the diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index 6967368e8d..3056a686c7 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -187,6 +187,11 @@ class MappingCartesian : public Mapping<dim,spacedim> */ Tensor<1,dim> length; + /** + * The volume element + */ + double volume_element; + /** * Vector of all quadrature * points. Especially, all diff --git a/deal.II/deal.II/source/fe/fe_abf.cc b/deal.II/deal.II/source/fe/fe_abf.cc index ce017586d3..55d99e23e3 100644 --- a/deal.II/deal.II/source/fe/fe_abf.cc +++ b/deal.II/deal.II/source/fe/fe_abf.cc @@ -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 @@ -56,7 +56,7 @@ FE_ABF<dim>::FE_ABF (const unsigned int deg) Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = this->contravariant; + this->mapping_type = mapping_piola; // First, initialize the // generalized support points and // quadrature weights, since they @@ -539,12 +539,12 @@ FE_ABF<dim>::update_each (const UpdateFlags flags) const if (flags & update_values) out |= update_values | update_covariant_transformation - | update_contravariant_transformation + | update_piola | update_cell_JxW_values | update_JxW_values; if (flags & update_gradients) out |= update_gradients | update_covariant_transformation - | update_contravariant_transformation + | update_piola | update_cell_JxW_values | update_JxW_values; //TODO: Set update flags appropriately and figure out, how the second 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 5b53f9613a..3ca14a3a4e 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -257,9 +257,10 @@ FE_PolyTensor<POLY,dim,spacedim>::shape_grad_grad_component (const unsigned int template <class POLY, int dim, int spacedim> typename Mapping<dim,spacedim>::InternalDataBase * -FE_PolyTensor<POLY,dim,spacedim>::get_data (const UpdateFlags update_flags, - const Mapping<dim,spacedim> &mapping, - const Quadrature<dim> &quadrature) const +FE_PolyTensor<POLY,dim,spacedim>::get_data ( + const UpdateFlags update_flags, + const Mapping<dim,spacedim> &mapping, + const Quadrature<dim> &quadrature) const { // generate a new data object and // initialize some fields @@ -371,11 +372,6 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( // possible InternalData &fe_data = dynamic_cast<InternalData &> (fedata); -// Assert(mapping_type == independent -// || ( mapping_type == independent_on_cartesian -// && dynamic_cast<const MappingCartesian<dim>*>(&mapping) != 0), -// ExcNotImplemented()); - const unsigned int n_q_points = quadrature.size(); const UpdateFlags flags(fe_data.current_update_flags()); @@ -400,8 +396,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( { switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) @@ -409,13 +404,13 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( break; } - case covariant: - case contravariant: + case mapping_covariant: + case mapping_contravariant: { // Use auxiliary vector for // transformation std::vector<Tensor<1,dim> > shape_values (n_q_points); - if (mapping_type == covariant) + if (mapping_type == mapping_covariant) mapping.transform(fe_data.shape_values[i], shape_values, mapping_data, mapping_covariant); else @@ -425,9 +420,12 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( // then copy over to target: 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 == contravariant ? + J = (mapping_type == mapping_contravariant ? data.JxW_values[k] / quadrature.weight(k) : 1.0); @@ -438,6 +436,18 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( } 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; + } default: Assert(false, ExcNotImplemented()); @@ -451,8 +461,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { mapping.transform(fe_data.shape_grads[i], shape_grads1, mapping_data, mapping_covariant); @@ -462,7 +471,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( break; } - case covariant: + case mapping_covariant: { mapping.transform(fe_data.shape_grads[i], shape_grads1, mapping_data, mapping_covariant); @@ -481,8 +490,8 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( break; } - - case contravariant: + case mapping_piola: + case mapping_contravariant: { mapping.transform(fe_data.shape_grads[i], shape_grads1, mapping_data, mapping_covariant); @@ -494,11 +503,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values ( for (unsigned int d=0; d<dim; ++d) { // Recompute determinant - const double - J = (mapping_type == contravariant ? - data.JxW_values[k] / quadrature.weight(k) - : - 1.0); + const double J = data.JxW_values[k] / quadrature.weight(k); data.shape_gradients[first+d][k] = sign_change[i] * shape_grads2[k][d] / J; } @@ -574,8 +579,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( { switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) @@ -583,13 +587,13 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( break; } - case covariant: - case contravariant: + case mapping_covariant: + case mapping_contravariant: { // Use auxiliary vector // for transformation std::vector<Tensor<1,dim> > shape_values (n_q_points); - if (mapping_type == covariant) + 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 @@ -611,7 +615,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( // face in // fill_fe_face_values const double - J = (mapping_type == contravariant ? + J = (mapping_type == mapping_contravariant ? data.cell_JxW_values[k] / quadrature.weight(k) : 1.0); @@ -623,6 +627,18 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( 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; + } default: Assert(false, ExcNotImplemented()); @@ -636,8 +652,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_covariant); @@ -647,7 +662,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( break; } - case covariant: + case mapping_covariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_covariant); @@ -665,8 +680,8 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( data.shape_gradients[first+d][k] = shape_grads1[k][d]; break; } - - case contravariant: + 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); @@ -680,8 +695,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values ( // recompute // determinant in the // same way as above - const double - J = data.cell_JxW_values[k] / quadrature.weight(k); + 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; } @@ -755,8 +769,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( { switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { for (unsigned int k=0; k<n_q_points; ++k) for (unsigned int d=0; d<dim; ++d) @@ -764,13 +777,13 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( break; } - case covariant: - case contravariant: + case mapping_covariant: + case mapping_contravariant: { // Use auxiliary vector for // transformation std::vector<Tensor<1,dim> > shape_values (n_q_points); - if (mapping_type == covariant) + 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 @@ -792,7 +805,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( // face in // fill_fe_face_values const double - J = (mapping_type == contravariant ? + J = (mapping_type == mapping_contravariant ? data.cell_JxW_values[k] / quadrature.weight(k) : 1.0); @@ -804,7 +817,18 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( 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; + } default: Assert(false, ExcNotImplemented()); } @@ -817,8 +841,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_covariant); @@ -828,7 +851,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( break; } - case covariant: + case mapping_covariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_covariant); @@ -846,8 +869,9 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( data.shape_gradients[first+d][k] = shape_grads1[k][d]; break; } - - case contravariant: + + 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); @@ -861,8 +885,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values ( // recompute // determinant in the // same way as above - const double - J = data.cell_JxW_values[k] / quadrature.weight(k); + 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; } @@ -913,8 +936,7 @@ template <class POLY, int dim, int spacedim> UpdateFlags FE_PolyTensor<POLY,dim,spacedim>::update_once (const UpdateFlags flags) const { - Assert (mapping_type != no_mapping, ExcNotInitialized()); - const bool values_once = (mapping_type == independent); + const bool values_once = (mapping_type == mapping_none); UpdateFlags out = update_default; @@ -929,8 +951,7 @@ template <class POLY, int dim, int spacedim> UpdateFlags FE_PolyTensor<POLY,dim,spacedim>::update_each (const UpdateFlags flags) const { - Assert (mapping_type != no_mapping, ExcNotInitialized()); - const bool values_once = (mapping_type == independent); + const bool values_once = (mapping_type == mapping_none); UpdateFlags out = update_default; 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 6b0b24da4e..900ebb1c15 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -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 @@ -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 = this->contravariant; + this->mapping_type = mapping_piola; // First, initialize the // generalized support points and // quadrature weights, since they @@ -464,12 +464,12 @@ FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const if (flags & update_values) out |= update_values | update_covariant_transformation - | update_contravariant_transformation + | update_piola | update_cell_JxW_values | update_JxW_values; if (flags & update_gradients) out |= update_gradients | update_covariant_transformation - | update_contravariant_transformation + | update_piola | update_cell_JxW_values | update_JxW_values; if (flags & update_hessians) 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 bf20b13275..8af15830ad 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 @@ -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 @@ -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 = this->independent_on_cartesian; + this->mapping_type = mapping_none; // These must be done first, since // they change the evaluation of // basis functions diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 9613be6852..c9f4dbb9c3 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -86,7 +86,7 @@ MappingCartesian<dim, spacedim>::get_data (const UpdateFlags update_flags, data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); data->update_flags = data->update_once | data->update_each; - + return data; } @@ -103,7 +103,7 @@ MappingCartesian<dim, spacedim>::get_face_data (const UpdateFlags update_flags, data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); data->update_flags = data->update_once | data->update_each; - + return data; } @@ -120,7 +120,7 @@ MappingCartesian<dim, spacedim>::get_subface_data (const UpdateFlags update_flag data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); data->update_flags = data->update_once | data->update_each; - + return data; } @@ -325,13 +325,15 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator& cell, // equal and are the product of the // local lengths in each coordinate // direction - if (data.current_update_flags() & update_JxW_values) + if (data.current_update_flags() & (update_JxW_values | update_volume_elements)) { double J = data.length[0]; for (unsigned int d=1;d<dim;++d) J *= data.length[d]; - for (unsigned int i=0; i<JxW_values.size();++i) - JxW_values[i] = J * q.weight(i); + data.volume_element = J; + if (data.current_update_flags() & update_JxW_values) + for (unsigned int i=0; i<JxW_values.size();++i) + JxW_values[i] = J * q.weight(i); } // "compute" Jacobian at the quadrature // points, which are all the same @@ -402,7 +404,13 @@ MappingCartesian<dim, spacedim>::fill_fe_face_values (const typename Triangulati for (unsigned int i=0; i<boundary_forms.size();++i) boundary_forms[i] = J * normal_vectors[i]; - //TODO: Implement cell_JxW_values + if (data.current_update_flags() & update_volume_elements) + { + J = data.length[0]; + for (unsigned int d=1;d<dim;++d) + J *= data.length[d]; + data.volume_element = J; + } } @@ -464,7 +472,13 @@ MappingCartesian<dim, spacedim>::fill_fe_subface_values (const typename Triangul for (unsigned int i=0; i<boundary_forms.size();++i) boundary_forms[i] = J * normal_vectors[i]; - //TODO: Implement cell_JxW_values + if (data.current_update_flags() & update_volume_elements) + { + J = data.length[0]; + for (unsigned int d=1;d<dim;++d) + J *= data.length[d]; + data.volume_element = J; + } } @@ -541,6 +555,19 @@ MappingCartesian<dim,spacedim>::transform ( 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; + } default: Assert(false, ExcNotImplemented()); } diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 5c0b52828c..d229bda802 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -393,6 +393,7 @@ MappingQ1<dim,spacedim>::update_each (const UpdateFlags in) const | update_cell_JxW_values | update_boundary_forms | update_normal_vectors + | update_volume_elements | update_jacobians | update_jacobian_grads | update_inverse_jacobians)); @@ -478,6 +479,9 @@ MappingQ1<dim,spacedim>::compute_data (const UpdateFlags update_flags, if (flags & update_contravariant_transformation) data.contravariant.resize(n_original_q_points); + if (flags & update_volume_elements) + data.volume_elements.resize(n_original_q_points); + if (flags & update_jacobian_grads) data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points); @@ -675,6 +679,9 @@ MappingQ1<dim,spacedim>::compute_fill (const typename Triangulation<dim,spacedim data.contravariant[point][i][j] += data_derv[j] * supp_pts[i]; } + if (update_flags & update_volume_elements) + for (unsigned int point=0; point<n_q_points; ++point) + data.volume_elements[point] = determinant(data.contravariant[point]); // invert contravariant for // covariant transformation @@ -1204,6 +1211,22 @@ MappingQ1<dim,spacedim>::transform ( } 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; + } default: Assert(false, ExcNotImplemented()); } -- 2.39.5