From f150c5c5f88bf5eafb1497064f7926fd01bf241e Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 18 Sep 2002 15:00:11 +0000 Subject: [PATCH] Add support for face and subface values. Implement a few other missing things. Fix some bugs. git-svn-id: https://svn.dealii.org/trunk@6452 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_nedelec.cc | 389 ++++++++++++++++----- deal.II/deal.II/source/fe/fe_nedelec_1d.cc | 4 +- deal.II/deal.II/source/fe/fe_nedelec_3d.cc | 2 - 3 files changed, 298 insertions(+), 97 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index 2f283743cf..2db8b0705f 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -874,9 +874,9 @@ FE_Nedelec::fill_fe_values (const Mapping &mapping, std::vector > shape_grads1 (n_q_points); std::vector > shape_grads2 (n_q_points); - Assert (data.shape_values.n_rows() == dofs_per_cell * dim, + Assert (data.shape_gradients.n_rows() == dofs_per_cell * dim, ExcInternalError()); - Assert (data.shape_values.n_cols() == n_q_points, + Assert (data.shape_gradients.n_cols() == n_q_points, ExcInternalError()); // loop over all shape @@ -934,107 +934,295 @@ FE_Nedelec::fill_fe_values (const Mapping &mapping, compute_2nd (mapping, cell, 0, mapping_data, fe_data, data); fe_data.first_cell = false; -} +}; template void -FE_Nedelec::fill_fe_face_values (const Mapping &/*mapping*/, - const typename DoFHandler::cell_iterator &/*cell*/, - const unsigned int /*face*/, - const Quadrature &/*quadrature*/, - typename Mapping::InternalDataBase &/*mapping_data*/, - typename Mapping::InternalDataBase &/*fedata*/, - FEValuesData &/*data*/) const +FE_Nedelec::fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { -//TODO!! -// // convert data object to internal -// // data for this class. fails with -// // an exception if that is not -// // possible -// InternalData &fe_data = dynamic_cast (fedata); - -// // offset determines which data set -// // to take (all data sets for all -// // faces are stored contiguously) -// const unsigned int offset = face * quadrature.n_quadrature_points; - -// const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &fe_data = dynamic_cast (fedata); + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + const unsigned int offset = face * quadrature.n_quadrature_points; + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + Assert (fe_data.shape_values.n_cols() == + GeometryInfo::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_values (n_q_points); -// for (unsigned int k=0; k::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; k void -FE_Nedelec::fill_fe_subface_values (const Mapping &/*mapping*/, - const typename DoFHandler::cell_iterator &/*cell*/, - const unsigned int /*face*/, - const unsigned int /*subface*/, - const Quadrature &/*quadrature*/, - typename Mapping::InternalDataBase &/*mapping_data*/, - typename Mapping::InternalDataBase &/*fedata*/, - FEValuesData &/*data*/) const +FE_Nedelec::fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const unsigned int subface, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { -//TODO!! -// // convert data object to internal -// // data for this class. fails with -// // an exception if that is not -// // possible -// InternalData &fe_data = dynamic_cast (fedata); - -// // offset determines which data set -// // to take (all data sets for all -// // sub-faces are stored contiguously) -// const unsigned int offset = (face * GeometryInfo::subfaces_per_face + subface) -// * quadrature.n_quadrature_points; - -// const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - -// for (unsigned int k=0; k (fedata); + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + const unsigned int offset = ((face * GeometryInfo::subfaces_per_face + subface) + * quadrature.n_quadrature_points); + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + Assert (fe_data.shape_values.n_cols() == + GeometryInfo::faces_per_cell * n_q_points, + ExcInternalError()); -// if (flags & update_gradients) -// { -// Assert (data.shape_gradients[k].size() + offset<= -// fe_data.shape_gradients[k].size(), -// ExcInternalError()); -// mapping.transform_covariant(data.shape_gradients[k], -// fe_data.shape_gradients[k], -// mapping_data, offset); -// }; -// } + std::vector > shape_values (n_q_points); + + Assert (data.shape_values.n_rows() == dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_values.n_cols() == n_q_points, + ExcInternalError()); + + for (unsigned int k=0; k::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; k::has_support_on_face (const unsigned int shape_index, ExcIndexRange (shape_index, 0, this->dofs_per_cell)); Assert (face_index < GeometryInfo::faces_per_cell, ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); - -//TODO: fix for higher orders. correct now for lowest order, all dimensions -//TODO!! -// can this be done in a way that is dimension and degree independent? + + switch (degree) + { + case 1: + { + // only on non-adjacent faces + // are the values actually + // zero. list these in a + // table + const unsigned int opposite_faces_2d[GeometryInfo<2>::faces_per_cell] + = { 2, 3, 0, 1}; + const unsigned int opposite_faces_3d[GeometryInfo<3>::faces_per_cell] + = { 1, 0, 4, 5, 2, 3}; + switch (dim) + { + case 2: return (face_index != opposite_faces_2d[shape_index]); + case 3: return (face_index != opposite_faces_3d[shape_index]); + default: Assert (false, ExcNotImplemented()); + }; + }; + + default: // other degree + Assert (false, ExcNotImplemented()); + }; - // all degrees of freedom are on - // lines, so also on a face. the - // question is whether it has - // support on this particular face - Assert (false, ExcNotImplemented()); return true; } diff --git a/deal.II/deal.II/source/fe/fe_nedelec_1d.cc b/deal.II/deal.II/source/fe/fe_nedelec_1d.cc index 10cf287383..1babaddbd8 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec_1d.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec_1d.cc @@ -13,8 +13,8 @@ -// only compile this file if in 1d. note that Nedelec elemets do not -// make much sense in 1d, so this file only contains dummy +// only compile this file if in 1d. note that Nedelec elements do not +// make much sense in 1d anyway, so this file only contains dummy // implementations to avoid linker errors due to missing symbols #if deal_II_dimension == 1 diff --git a/deal.II/deal.II/source/fe/fe_nedelec_3d.cc b/deal.II/deal.II/source/fe/fe_nedelec_3d.cc index 38638df49f..6afa6fbe65 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec_3d.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec_3d.cc @@ -24,8 +24,6 @@ // the child cell are obtained from the degrees of freedom on the // mother cell // -// TODO: [Anna] check whether the following paragraph is correct. if so, then please multiply the values in the eight following matrices by two - // note the following: since the shape functions themselves and not // only the gradients are transformed using the mapping object from // the unit cell to the real cell, the actual values of the function -- 2.39.5