From f936cc2d776bd73db953aa43458dec339cb856ab Mon Sep 17 00:00:00 2001 From: kayser-herold Date: Tue, 22 Aug 2006 17:39:31 +0000 Subject: [PATCH] Extended the mapping classes with functionality to provide the determinant of the Jacobian of the cell mapping. In particular this functionality is required for the fe_poly_tensor class, where a Piola transformation is used. Furthermore the procedure to determine required changes of the sign of the Face DoFs of RT and ABF elements was extended to be able to cope with refined grids. The last patch implements the third (difficult) case for computing the interface contraints in the hp case. git-svn-id: https://svn.dealii.org/trunk@13717 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_poly_tensor.h | 16 +++ deal.II/deal.II/include/fe/fe_update_flags.h | 14 +- deal.II/deal.II/include/fe/fe_values.h | 25 ++++ deal.II/deal.II/include/fe/mapping.h | 6 +- .../deal.II/include/fe/mapping_cartesian.h | 12 +- deal.II/deal.II/include/fe/mapping_q.h | 6 +- deal.II/deal.II/include/fe/mapping_q1.h | 18 ++- deal.II/deal.II/include/grid/tria_accessor.h | 5 + deal.II/deal.II/source/dofs/dof_tools.cc | 102 +++++++++++++- deal.II/deal.II/source/fe/fe_abf.cc | 5 +- deal.II/deal.II/source/fe/fe_poly_tensor.cc | 130 +++++++++++------- .../deal.II/source/fe/fe_raviart_thomas.cc | 2 + deal.II/deal.II/source/fe/fe_system.cc | 1 + deal.II/deal.II/source/fe/fe_values.cc | 10 +- .../deal.II/source/fe/mapping_cartesian.cc | 16 ++- deal.II/deal.II/source/fe/mapping_q.cc | 16 ++- deal.II/deal.II/source/fe/mapping_q1.cc | 51 +++++-- 17 files changed, 341 insertions(+), 94 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 efd05294a1..577cd571e8 100644 --- a/deal.II/deal.II/include/fe/fe_poly_tensor.h +++ b/deal.II/deal.II/include/fe/fe_poly_tensor.h @@ -261,6 +261,22 @@ class FE_PolyTensor : public FiniteElement contravariant }; + + + /** + * On noncartesian grids, the + * sign of the DoFs associated with + * the faces of the elements has + * to be changed in some cases. + * This procedure implements an + * algorithm, which determines + * the DoFs, which need this + * sign change for a given cell. + */ + void get_face_sign_change (const typename Triangulation::cell_iterator &cell, + std::vector &face_sign) const; + + /** * The mapping type to be used to * map shape functions from the 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 5bec6b407f..a83e7634ea 100644 --- a/deal.II/deal.II/include/fe/fe_update_flags.h +++ b/deal.II/deal.II/include/fe/fe_update_flags.h @@ -199,7 +199,19 @@ enum UpdateFlags * transformation defined by * the Mapping. */ - update_transformation_gradients = 0x1000 + update_transformation_gradients = 0x1000, + + /** + * Compute the JxW values + * on faces for the cell mapping + * and not the face mapping! + * This functionality is required + * for the Piola mapping which + * is used in conjunction with + * H_div subspaces like RT and ABF. + */ + update_cell_JxW_values = 0x2000 + }; diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 0b6e2d5a7f..5941ecd692 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -183,6 +183,31 @@ class FEValuesData */ std::vector JxW_values; + /** + * Store an array of weights + * times the Jacobi determinant + * at the quadrature points. This + * function is reset each time + * reinit() is called. The + * Jacobi determinant is actually + * the reciprocal value of the + * Jacobi matrices stored in this + * class, see the general + * documentation of this class + * for more information. + * + * In contrast to JxW_values, + * this vector uses the transformation + * of the cell even in case the + * corresponding FEValues object + * refers to a face. Currently the + * only purpose is to provide + * exactly this value to the + * transformation which is used + * for the PolyTensor elements. + */ + std::vector cell_JxW_values; + /** * Array of quadrature points. This array * is set up upon calling reinit() and diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index e4bfa3aaca..bb0660c0e3 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -472,7 +472,8 @@ class Mapping : public Subscriptor std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_form, - std::vector > &normal_vectors) const = 0; + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const = 0; /** * See above. @@ -486,7 +487,8 @@ class Mapping : public Subscriptor std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_form, - std::vector > &normal_vectors) const = 0; + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const = 0; /** * Give class @p FEValues access diff --git a/deal.II/deal.II/include/fe/mapping_cartesian.h b/deal.II/deal.II/include/fe/mapping_cartesian.h index 804ab76a9e..967b58c862 100644 --- a/deal.II/deal.II/include/fe/mapping_cartesian.h +++ b/deal.II/deal.II/include/fe/mapping_cartesian.h @@ -68,7 +68,8 @@ class MappingCartesian : public Mapping std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_form, - std::vector > &normal_vectors) const ; + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const ; virtual void fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, @@ -78,7 +79,8 @@ class MappingCartesian : public Mapping std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_form, - std::vector > &normal_vectors) const ; + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const ; /** * Implementation of the @@ -232,7 +234,8 @@ template <> void MappingCartesian<1>::fill_fe_face_values ( std::vector >&, std::vector&, std::vector >&, - std::vector >&) const; + std::vector >&, + std::vector&) const; template <> void MappingCartesian<1>::fill_fe_subface_values ( const Triangulation<1>::cell_iterator &, @@ -243,7 +246,8 @@ template <> void MappingCartesian<1>::fill_fe_subface_values ( std::vector >&, std::vector&, std::vector >&, - std::vector >&) const; + std::vector >&, + std::vector&) const; #endif // DOXYGEN diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index f6f4d01ba8..a98a841d53 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -220,7 +220,8 @@ class MappingQ : public MappingQ1 typename std::vector > &quadrature_points, std::vector &JxW_values, typename std::vector > &exterior_form, - typename std::vector > &normal_vectors) const ; + typename std::vector > &normal_vectors, + std::vector &cell_JxW_values) const ; /** * Implementation of the interface in @@ -235,7 +236,8 @@ class MappingQ : public MappingQ1 typename std::vector > &quadrature_points, std::vector &JxW_values, typename std::vector > &exterior_form, - typename std::vector > &normal_vectors) const ; + typename std::vector > &normal_vectors, + std::vector &cell_JxW_values) const ; /** * For dim=2,3. Append the diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index 2a9fa19e08..415e236fa4 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -397,7 +397,8 @@ class MappingQ1 : public Mapping typename std::vector > &quadrature_points, std::vector &JxW_values, typename std::vector > &boundary_form, - typename std::vector > &normal_vectors) const ; + typename std::vector > &normal_vectors, + std::vector &cell_JxW_values) const ; /** * Implementation of the interface in @@ -412,7 +413,8 @@ class MappingQ1 : public Mapping typename std::vector > &quadrature_points, std::vector &JxW_values, typename std::vector > &boundary_form, - typename std::vector > &normal_vectors) const ; + typename std::vector > &normal_vectors, + std::vector &cell_JxW_values) const ; /** * Compute shape values and/or @@ -483,7 +485,8 @@ class MappingQ1 : public Mapping std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_form, - std::vector > &normal_vectors) const; + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const; /** * Compute shape values and/or @@ -712,7 +715,8 @@ template <> void MappingQ1<1>::compute_fill_face ( std::vector > &, std::vector &, std::vector > &, - std::vector > &) const; + std::vector > &, + std::vector &) const; template <> void MappingQ1<1>::fill_fe_face_values ( const Triangulation<1>::cell_iterator &, @@ -722,7 +726,8 @@ template <> void MappingQ1<1>::fill_fe_face_values ( std::vector >&, std::vector&, std::vector >&, - std::vector >&) const; + std::vector >&, + std::vector&) const; template <> void MappingQ1<1>::fill_fe_subface_values ( const Triangulation<1>::cell_iterator &, @@ -733,7 +738,8 @@ template <> void MappingQ1<1>::fill_fe_subface_values ( std::vector >&, std::vector&, std::vector >&, - std::vector >&) const; + std::vector >&, + std::vector&) const; #endif // DOXYGEN diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index 06a8834568..de217b11d8 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -2538,6 +2538,11 @@ class CellAccessor : public TriaObjectAccessor * @p ith neighbor. If the * neighbor does not exist, an * invalid iterator is returned. + * + * Note (cf. TriaLevel<0>): + * The neighbor of a cell has at most the + * same level as this cell, i.e. it may + * or may not be refined. */ TriaIterator > neighbor (const unsigned int i) const; diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 1e3b832cf7..d7a231fb65 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -2287,13 +2287,29 @@ namespace internal // that the other fine elements will // be constrained to the coarse element. { - Assert (false, ExcNotImplemented ()); + // Assert (false, ExcNotImplemented ()); typename DH::active_cell_iterator neighbor_child = cell->neighbor_child_on_subface (face, min_degree_subface); const unsigned int n_dofs_on_children = neighbor_child->get_fe().dofs_per_face; dofs_on_children.resize (n_dofs_on_children); + + // Get DoFs on child cell with + // lowest polynomial degree. + // All other DoFs will be constrained + // to the DoFs of this face. + const unsigned int subface_fe_index + = neighbor_child->active_fe_index(); + + // Same procedure as for the + // mother cell. Extract the face + // DoFs from the cell DoFs. + cell->face(face)->child(min_degree_subface) + ->get_dof_indices (dofs_on_children, + subface_fe_index); + + // The idea is to introduce // a "virtual" intermediate coarse // level face with the lowest @@ -2361,12 +2377,88 @@ namespace internal // Now create the element // constraint for this subface. - FullMatrix face_constraints_m (n_dofs_on_children, - n_dofs_on_children); + FullMatrix fc_sface_ipol (n_dofs_on_children, + n_dofs_on_children); + FullMatrix fc_ipol_sface (n_dofs_on_children, + n_dofs_on_children); neighbor_child->get_fe().get_subface_interpolation_matrix (neighbor_child->get_fe (), min_degree_subface, - face_constraints_m); -//TODO: Continue ... + fc_sface_ipol); + // Invert it, to get a mapping from the DoFs of the + // "Master-subface" to the intermediate layer. + fc_ipol_sface.invert (fc_sface_ipol); + + // Create constraint matrix for the mother face. + const unsigned int n_dofs_on_mother = cell->get_fe().dofs_per_face; + dofs_on_mother.resize (n_dofs_on_mother); + cell->face(face)->get_dof_indices (dofs_on_mother, cell->active_fe_index ()); + + FullMatrix fc_mother_ipol (n_dofs_on_children, + n_dofs_on_mother); + FullMatrix fc_mother_sface (n_dofs_on_children, + n_dofs_on_mother); + neighbor_child->get_fe ().get_face_interpolation_matrix (cell->get_fe(), + fc_mother_ipol); + fc_ipol_sface.mmult (fc_mother_sface, fc_mother_ipol); + + // Add constraints to global constraint + // matrix. + filter_constraints (dofs_on_children, + dofs_on_mother, + fc_mother_sface, + constraints); + + // Now create constraint matrices for + // the subfaces and assemble them + for (unsigned int c=0; c::subfaces_per_face; ++c) + { + // As the "Master-subface" does not need constraints, skip it. + if (c != min_degree_subface) + { + typename DH::active_cell_iterator neighbor_child_slave + = cell->neighbor_child_on_subface (face, c); + const unsigned int n_dofs_on_mother = neighbor_child_slave->get_fe().dofs_per_face; + dofs_on_mother.resize (n_dofs_on_mother); + + // Find face number on the finer + // neighboring cell, which is + // shared the face with the + // face of the coarser cell. + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face); + Assert (neighbor_child_slave->face(neighbor2) == cell->face(face)->child(c), + ExcInternalError()); + + // Same procedure as for the + // mother cell. Extract the face + // DoFs from the cell DoFs. + const unsigned int subface_fe_index + = neighbor_child_slave->active_fe_index(); + + cell->face(face)->child(c) + ->get_dof_indices (dofs_on_mother, + subface_fe_index); + + // Now create the element + // constraint for this subface. + FullMatrix fc_child_sface_ipol (n_dofs_on_children, + n_dofs_on_mother); + FullMatrix fc_child_sface_sface (n_dofs_on_children, + n_dofs_on_mother); + neighbor_child->get_fe ().get_subface_interpolation_matrix + (neighbor_child_slave->get_fe(), + c, fc_child_sface_ipol); + + fc_ipol_sface.mmult (fc_child_sface_sface, fc_child_sface_ipol); + + // Add constraints to global constraint + // matrix. + filter_constraints (dofs_on_children, + dofs_on_mother, + fc_child_sface_sface, + constraints); + } + } } } else diff --git a/deal.II/deal.II/source/fe/fe_abf.cc b/deal.II/deal.II/source/fe/fe_abf.cc index c16dc5a6fa..ea3a692917 100644 --- a/deal.II/deal.II/source/fe/fe_abf.cc +++ b/deal.II/deal.II/source/fe/fe_abf.cc @@ -62,9 +62,6 @@ FE_ABF::FE_ABF (const unsigned int deg) FullMatrix M(n_dofs, n_dofs); FETools::compute_node_matrix(M, *this); - //TODO: Remove debugging output - // M.print (std::cout); - this->inverse_node_matrix.reinit(n_dofs, n_dofs); this->inverse_node_matrix.invert(M); // From now on, the shape functions @@ -532,10 +529,12 @@ FE_ABF::update_each (const UpdateFlags flags) const if (flags & update_values) out |= update_values | update_covariant_transformation | update_contravariant_transformation + | update_cell_JxW_values | update_JxW_values; if (flags & update_gradients) out |= update_gradients | update_covariant_transformation | update_contravariant_transformation + | update_cell_JxW_values | update_JxW_values; //TODO: Set update flags appropriately and figure out, how the second // derivatives for the RT elements can be computed correctly. 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 fa8812ced2..cf97b95d1e 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -154,6 +154,76 @@ FE_PolyTensor::shape_grad_grad_component ( +//--------------------------------------------------------------------------- +// Utility method, which is used to determine the change of sign for +// the DoFs on the faces of the given cell. +//--------------------------------------------------------------------------- + +template +void +FE_PolyTensor::get_face_sign_change ( + const typename Triangulation::cell_iterator &cell, + std::vector &face_sign) const +{ + Assert (face_sign.size () == this->dofs_per_cell, + ExcDimensionMismatch (face_sign.size (), this->dofs_per_cell)); + + // Default is no sign change. I.e. multiply by one. + std::fill (face_sign.begin (), face_sign.end (), 1.0); + +#if deal_II_dimension > 1 + const unsigned int dofs_per_face = this->dofs_per_face; + + if (dim == 2) + { + for (unsigned int f = GeometryInfo::faces_per_cell / 2; + f < GeometryInfo::faces_per_cell; ++f) + { + typename Triangulation::face_iterator face = cell->face (f); + if (!face->at_boundary ()) + { + const unsigned int neighbor_level = cell->neighbor (f)->level (); + const unsigned int cell_level = cell->level (); + + unsigned int nn = (unsigned int) -1; + + // Same level, the easy case. + if (neighbor_level == cell_level) + nn = cell->neighbor_of_neighbor (f); + else + { + // Neighbor is more refined (due to the internal restrictions, + // this case should never be encountered! + if (neighbor_level > cell_level) + { + Assert (false, + ExcMessage ("neighbor_level is larger than cell level!")); + } + // Neighbor is less refined + else + nn = cell->neighbor_of_coarser_neighbor (f).first; + } + + Assert (nn != (unsigned int) -1, + ExcInternalError ()); + + if (nn < GeometryInfo::faces_per_cell / 2) + { + for (unsigned int j = 0; j < dofs_per_face; ++j) + face_sign[f * dofs_per_face + j] = -1.0; + } + } + } + } + else + { + // TODO: Think about 3D!. + } +#endif + +} + + //--------------------------------------------------------------------------- // Data field initialization //--------------------------------------------------------------------------- @@ -289,32 +359,7 @@ FE_PolyTensor::fill_fe_values ( // Compute eventual sign changes depending on the neighborhood // between two faces. std::vector sign_change (n_dofs, 1.0); - -#if deal_II_dimension > 1 - const unsigned int dofs_per_face = this->dofs_per_face; - - if (dim == 2) - { - for (unsigned int f = GeometryInfo::faces_per_cell / 2; - f < GeometryInfo::faces_per_cell; ++f) - { - typename Triangulation::face_iterator face = cell->face (f); - if (!face->at_boundary ()) - { - unsigned int nn = cell->neighbor_of_neighbor (f); - if (nn < GeometryInfo::faces_per_cell / 2) - { - for (unsigned int j = 0; j < dofs_per_face; ++j) - sign_change[f * dofs_per_face + j] = -1.0; - } - } - } - } - else - { - // TODO: Think about 3D!. - } -#endif + get_face_sign_change (cell, sign_change); for (unsigned int i=0; i::fill_fe_face_values ( // Compute eventual sign changes depending on the neighborhood // between two faces. std::vector sign_change (n_dofs, 1.0); -#if deal_II_dimension > 1 - const unsigned int dofs_per_face = this->dofs_per_face; - - if (dim == 2) - { - for (unsigned int f = GeometryInfo::faces_per_cell / 2; - f < GeometryInfo::faces_per_cell; ++f) - { - typename Triangulation::face_iterator face = cell->face (f); - if (!face->at_boundary ()) - { - unsigned int nn = cell->neighbor_of_neighbor (f); - if (nn < GeometryInfo::faces_per_cell / 2) - { - for (unsigned int j = 0; j < dofs_per_face; ++j) - sign_change[f * dofs_per_face + j] = -1.0; - } - } - } - } - else - { - // TODO: Think about 3D!. - } -#endif + get_face_sign_change (cell, sign_change); for (unsigned int i=0; i::fill_fe_face_values ( //TODO: Note that the Jacobian we compute here is that of the transformation from the reference FACE to the real FACE. However, what we need is the Jacobian of the transformation of the reference to the real CELL, which differs by one power of h. Consequently, the tests fe/rt_{11,13,14,15} presently compute the wrong results double J = 1.0; if (mapping_type == contravariant) - J = data.JxW_values[k] / quadrature.weight(k); + J = data.cell_JxW_values[k] / quadrature.weight(k); for (unsigned int d=0; d::fill_fe_face_values ( { // Recompute determinant //TODO: Note that the Jacobian we compute here is that of the transformation from the reference FACE to the real FACE. However, what we need is the Jacobian of the transformation of the reference to the real CELL, which differs by one power of h. Consequently, the tests fe/rt_{11,13,14,15} presently compute the wrong results - double J = data.JxW_values[k] / quadrature.weight(k); + + double J = data.cell_JxW_values[k] / quadrature.weight(k); data.shape_gradients[first+d][k] = sign_change[i] * shape_grads2[k][d] / J; } @@ -639,6 +661,8 @@ FE_PolyTensor::fill_fe_subface_values ( // && dynamic_cast*>(&mapping) != 0), // ExcNotImplemented()); //TODO: Size assertions + +//TODO: Sign change for the face DoFs! for (unsigned int i=0; i::update_each (const UpdateFlags flags) const if (flags & update_values) out |= update_values | update_covariant_transformation | update_contravariant_transformation + | update_cell_JxW_values | update_JxW_values; if (flags & update_gradients) out |= update_gradients | update_covariant_transformation | update_contravariant_transformation + | update_cell_JxW_values | update_JxW_values; //TODO: Set update flags appropriately and figure out, how the second // derivatives for the RT elements can be computed correctly. diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index cd4f6982d3..0509aad8ec 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -832,6 +832,7 @@ compute_fill (const Mapping &mapping, // the reference element to the global // coordinate system. base_data.JxW_values = data.JxW_values; + base_data.cell_JxW_values = data.cell_JxW_values; // Make sure that in the diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index ad6160d4a9..9f7e38aaf1 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -312,6 +312,9 @@ FEValuesData::initialize (const unsigned int n_quadrature_points, if (flags & update_normal_vectors) this->normal_vectors.resize(n_quadrature_points); + + if (flags & update_cell_JxW_values) + this->cell_JxW_values.resize(n_quadrature_points); } @@ -1011,6 +1014,7 @@ FEValuesBase::memory_consumption () const MemoryConsumption::memory_consumption (this->quadrature_points) + MemoryConsumption::memory_consumption (this->normal_vectors) + MemoryConsumption::memory_consumption (this->boundary_forms) + + MemoryConsumption::memory_consumption (this->cell_JxW_values) + sizeof(this->update_flags) + MemoryConsumption::memory_consumption (n_quadrature_points) + MemoryConsumption::memory_consumption (dofs_per_cell) + @@ -1520,7 +1524,8 @@ void FEFaceValues::do_reinit (const unsigned int face_no) this->quadrature_points, this->JxW_values, this->boundary_forms, - this->normal_vectors); + this->normal_vectors, + this->cell_JxW_values); this->get_fe().fill_fe_face_values(this->get_mapping(), *this->present_cell, face_no, @@ -1770,7 +1775,8 @@ void FESubfaceValues::do_reinit (const unsigned int face_no, this->quadrature_points, this->JxW_values, this->boundary_forms, - this->normal_vectors); + this->normal_vectors, + this->cell_JxW_values); this->get_fe().fill_fe_subface_values(this->get_mapping(), *this->present_cell, diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 35e43ed841..58a0d82229 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -322,7 +322,8 @@ MappingCartesian::fill_fe_face_values (const typename Triangulation::c std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors, + std::vector &/*cell_JxW_values*/) const { // convert data object to internal // data for this class. fails with @@ -351,6 +352,8 @@ MappingCartesian::fill_fe_face_values (const typename Triangulation::c if (data.current_update_flags() & update_boundary_forms) for (unsigned int i=0; i::fill_fe_subface_values (const typename Triangulation std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors, + std::vector &/*cell_JxW_values*/) const { // convert data object to internal // data for this class. fails with @@ -394,6 +398,8 @@ MappingCartesian::fill_fe_subface_values (const typename Triangulation if (data.current_update_flags() & update_boundary_forms) for (unsigned int i=0; i::fill_fe_face_values (const Triangulation<1>::cell_iterator std::vector >&, std::vector&, std::vector >&, - std::vector >&) const + std::vector >&, + std::vector &/*cell_JxW_values*/) const { Assert(false, ExcNotImplemented()); } @@ -424,7 +431,8 @@ MappingCartesian<1>::fill_fe_subface_values (const Triangulation<1>::cell_iterat std::vector >&, std::vector&, std::vector >&, - std::vector >&) const + std::vector >&, + std::vector &/*cell_JxW_values*/) const { 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 c715dc1376..794c127720 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -342,7 +342,8 @@ MappingQ::fill_fe_face_values (const typename Triangulation::cell_iter std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &exterior_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const { // convert data object to internal // data for this class. fails with @@ -385,7 +386,10 @@ MappingQ::fill_fe_face_values (const typename Triangulation::cell_iter q.get_weights(), *p_data, quadrature_points, JxW_values, - exterior_forms, normal_vectors); + exterior_forms, normal_vectors, + cell_JxW_values); + + // TODO: Verify implementation of cell_JxW_values } @@ -399,7 +403,8 @@ MappingQ::fill_fe_subface_values (const typename Triangulation::cell_i std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &exterior_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const { // convert data object to internal // data for this class. fails with @@ -443,7 +448,10 @@ MappingQ::fill_fe_subface_values (const typename Triangulation::cell_i q.get_weights(), *p_data, quadrature_points, JxW_values, - exterior_forms, normal_vectors); + exterior_forms, normal_vectors, + cell_JxW_values); + + // TODO: Verify implementation of cell_JxW_values ... } diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index e4acb4ce35..b5bafeccbf 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -249,6 +249,7 @@ MappingQ1::update_each (const UpdateFlags in) const | update_covariant_transformation | update_contravariant_transformation | update_JxW_values + | update_cell_JxW_values | update_boundary_forms | update_normal_vectors)); @@ -271,6 +272,15 @@ MappingQ1::update_each (const UpdateFlags in) const | update_normal_vectors)) out |= update_contravariant_transformation; + // The contravariant transformation + // is a Piola transformation, which + // requires the determinant of the + // Jacobi matrix of the transformation. + // Therefore these values have to + // updated for each cell. + if (out & update_contravariant_transformation) + out |= update_JxW_values | update_cell_JxW_values; + return out; } @@ -582,7 +592,8 @@ MappingQ1::compute_fill_face (const typename Triangulation::cell_itera std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const { compute_fill (cell, n_q_points, data_set, data, quadrature_points); @@ -658,6 +669,23 @@ MappingQ1::compute_fill_face (const typename Triangulation::cell_itera normal_vectors[i] = boundary_forms[i] / f; } } + + // If the Piola transformation is to be used, the + // new flag, update_cell_JxW_values has to be set. + // It triggers the creation of values for the determinant + // of the mapping function on + if (update_flags & update_cell_JxW_values) + { + Assert (cell_JxW_values.size() == n_q_points, + ExcDimensionMismatch(cell_JxW_values.size(), n_q_points)); + Assert (data.contravariant.size() == n_q_points, + ExcDimensionMismatch(data.contravariant.size(), n_q_points)); + Assert (weights.size() == n_q_points, + ExcDimensionMismatch(weights.size(), n_q_points)); + for (unsigned int point=0; point::fill_fe_face_values (const typename Triangulation::cell_ite std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const { InternalData *data_ptr = dynamic_cast (&mapping_data); Assert(data_ptr!=0, ExcInternalError()); @@ -688,7 +717,8 @@ MappingQ1::fill_fe_face_values (const typename Triangulation::cell_ite quadrature_points, JxW_values, boundary_forms, - normal_vectors); + normal_vectors, + cell_JxW_values); } @@ -702,7 +732,8 @@ MappingQ1::fill_fe_subface_values (const typename Triangulation::cell_ std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &boundary_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors, + std::vector &cell_JxW_values) const { InternalData *data_ptr = dynamic_cast (&mapping_data); Assert(data_ptr!=0, ExcInternalError()); @@ -720,7 +751,8 @@ MappingQ1::fill_fe_subface_values (const typename Triangulation::cell_ quadrature_points, JxW_values, boundary_forms, - normal_vectors); + normal_vectors, + cell_JxW_values); } @@ -738,7 +770,8 @@ MappingQ1<1>::compute_fill_face (const Triangulation<1>::cell_iterator &, std::vector > &, std::vector &, std::vector > &, - std::vector > &) const + std::vector > &, + std::vector &/*cell_JxW_values*/) const { Assert(false, ExcNotImplemented()); } @@ -753,7 +786,8 @@ MappingQ1<1>::fill_fe_face_values (const Triangulation<1>::cell_iterator &, std::vector >&, std::vector&, std::vector >&, - std::vector >&) const + std::vector >&, + std::vector &/*cell_JxW_values*/) const { Assert(false, ExcNotImplemented()); } @@ -769,7 +803,8 @@ MappingQ1<1>::fill_fe_subface_values (const Triangulation<1>::cell_iterator &, std::vector >&, std::vector&, std::vector >&, - std::vector >&) const + std::vector >&, + std::vector &/*cell_JxW_values*/) const { Assert(false, ExcNotImplemented()); } -- 2.39.5