From: kayser-herold Date: Sat, 29 Jul 2006 20:40:45 +0000 (+0000) Subject: Added functionality which corrects some errors in the X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=52862969a4a4bd989b368b613e7a8395d8a6c9d2;p=dealii-svn.git Added functionality which corrects some errors in the code for the Raviart-Thomas and ABF elements. The first correction is a sign change of the Face-DoFs, which is required to work on more complicated domains. This fix is currently only applied in the 2D case. The second correction is related to the coordinate transformation from the reference element to the global coordinate system. This transformation now becomes the Piola transformation, which according to Brezzi/Fortin 1992 is the correct one for subspaces of H_div. git-svn-id: https://svn.dealii.org/trunk@13504 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 efe91afaaf..6df6043420 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -282,6 +282,39 @@ FE_PolyTensor::fill_fe_values ( ExcDimensionMismatch(fe_data.shape_values.size(), n_dofs)); Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_quad, ExcDimensionMismatch(fe_data.shape_values[0].size(), n_quad)); + + // Create table with sign changes, due to the special structure of the RT elements. + // TODO: Preliminary hack to demonstrate the overall prinicple! + + // 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 for (unsigned int i=0; i::fill_fe_values ( // then copy over to target: for (unsigned int k=0; k::fill_fe_values ( for (unsigned int k=0; k::fill_fe_face_values ( // && dynamic_cast*>(&mapping) != 0), // ExcNotImplemented()); //TODO: Size assertions + + // Create table with sign changes, due to the special structure of the RT elements. + // TODO: Preliminary hack to demonstrate the overall prinicple! + + // 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 for (unsigned int i=0; i::fill_fe_face_values ( // then copy over to target: for (unsigned int k=0; k::fill_fe_face_values ( for (unsigned int d=0;d::FE_RaviartThomas (const unsigned int deg) Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = this->covariant; + this->mapping_type = this->contravariant; // First, initialize the // generalized support points and // quadrature weights, since they @@ -449,11 +449,17 @@ FE_RaviartThomas::update_each (const UpdateFlags flags) const UpdateFlags out = update_default; if (flags & update_values) - out |= update_values | update_covariant_transformation; + out |= update_values | update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values; if (flags & update_gradients) - out |= update_gradients | update_covariant_transformation; + out |= update_gradients | update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values; + //TODO: Set update flags appropriately and figure out, how the second + // derivatives for the RT elements can be computed correctly. if (flags & update_second_derivatives) - out |= update_second_derivatives | update_covariant_transformation; + out |= update_second_derivatives | update_contravariant_transformation; return out; } diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 0721506226..9a85175dd1 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -821,6 +821,19 @@ compute_fill (const Mapping &mapping, FEValuesData & base_data = fe_data.get_fe_values_data(base_no); + //TODO: Think about a smarter alternative + // Copy quadrature points. These + // are required for computing the + // determinant in the FEPolyTensor + // class. The determinant is + // one ingredient of the Piola + // transformation, which is applied + // to correctly map the RT space from + // the reference element to the global + // coordinate system. + base_data.JxW_values = data.JxW_values; + + // Make sure that in the // case of fill_fe_values // the data is only copied