From ec2319d145e33a813081b75126eb44be611a4a81 Mon Sep 17 00:00:00 2001 From: kanschat 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 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 class FEValues; template class FEFaceValues; template 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 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 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 */ 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::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::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::shape_grad_grad_component (const unsigned int template typename Mapping::InternalDataBase * -FE_PolyTensor::get_data (const UpdateFlags update_flags, - const Mapping &mapping, - const Quadrature &quadrature) const +FE_PolyTensor::get_data ( + const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature) const { // generate a new data object and // initialize some fields @@ -371,11 +372,6 @@ FE_PolyTensor::fill_fe_values ( // possible InternalData &fe_data = dynamic_cast (fedata); -// Assert(mapping_type == independent -// || ( mapping_type == independent_on_cartesian -// && dynamic_cast*>(&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::fill_fe_values ( { switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { for (unsigned int k=0; k::fill_fe_values ( break; } - case covariant: - case contravariant: + case mapping_covariant: + case mapping_contravariant: { // Use auxiliary vector for // transformation std::vector > 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::fill_fe_values ( // then copy over to target: for (unsigned int k=0; k::fill_fe_values ( } break; } + case mapping_piola: + if (true) + { + std::vector > shape_values (n_q_points); + mapping.transform(fe_data.shape_values[i], shape_values, + mapping_data, mapping_piola); + for (unsigned int k=0; k::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::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::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::fill_fe_values ( for (unsigned int d=0; d::fill_fe_face_values ( { switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { for (unsigned int k=0; k::fill_fe_face_values ( break; } - case covariant: - case contravariant: + case mapping_covariant: + case mapping_contravariant: { // Use auxiliary vector // for transformation std::vector > 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::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::fill_fe_face_values ( break; } + case mapping_piola: + if (true) + { + std::vector > 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::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::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::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::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::fill_fe_subface_values ( { switch (mapping_type) { - case independent: - case independent_on_cartesian: + case mapping_none: { for (unsigned int k=0; k::fill_fe_subface_values ( break; } - case covariant: - case contravariant: + case mapping_covariant: + case mapping_contravariant: { // Use auxiliary vector for // transformation std::vector > 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::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::fill_fe_subface_values ( break; } - + case mapping_piola: + if (true) + { + std::vector > 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::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::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::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::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 UpdateFlags FE_PolyTensor::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 UpdateFlags FE_PolyTensor::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::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::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::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::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::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::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::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::fill_fe_face_values (const typename Triangulati for (unsigned int i=0; i::fill_fe_subface_values (const typename Triangul for (unsigned int i=0; i::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::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i::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::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::compute_fill (const typename Triangulation::transform ( } return; } + case mapping_piola: + if (true) + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField()); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField()); + + for (unsigned int i=0; i