From: Guido Kanschat Date: Wed, 25 May 2005 09:14:10 +0000 (+0000) Subject: debugging RaviartThomasNodal X-Git-Tag: v8.0.0~13816 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=058f76e534e0e47c622debc64fa2d5f091e4169a;p=dealii.git debugging RaviartThomasNodal git-svn-id: https://svn.dealii.org/trunk@10740 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 d18bb445ed..fa50303cfc 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -228,10 +228,10 @@ FE_PolyTensor::get_data (const UpdateFlags update_flags, if (inverse_node_matrix.n_cols() == 0) for (unsigned int i=0; idofs_per_cell; ++i) data->shape_values[i][k] = values[i]; - else - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - data->shape_values[i][k] = inverse_node_matrix(j,i) * values[j]; + else + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + data->shape_values[i][k] += inverse_node_matrix(j,i) * values[j]; if (flags & update_gradients) if (inverse_node_matrix.n_cols() == 0) @@ -240,7 +240,7 @@ FE_PolyTensor::get_data (const UpdateFlags update_flags, else for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; jdofs_per_cell; ++j) - data->shape_grads[i][k] = inverse_node_matrix(j,i) * grads[j]; + data->shape_grads[i][k] += inverse_node_matrix(j,i) * grads[j]; } return data; } @@ -282,7 +282,7 @@ FE_PolyTensor::fill_fe_values ( for (unsigned int i=0; i& fe) MappingCartesian mapping; Quadrature support_quadrature(points, weights); - UpdateFlags flags = update_values; + UpdateFlags flags = update_values | update_gradients; FEValues vals(mapping, fe, support_quadrature, flags); vals.reinit(dof.begin_active()); @@ -62,13 +62,15 @@ check_support_points (const FiniteElement& fe) for (unsigned int d=0;d 1.e-12) -// deallog << "Error" << vals.shape_value_component(i,k,d) << std::endl; + const double diff = std::abs(sf - vals.shape_value_component(i,k,d)); + if (diff > 1.e-12) + deallog << "Error values" << std::endl; + Tensor<1,dim> grad = fe.shape_grad_component(i,p,d); + grad -= vals.shape_grad_component(i,k,d); + if (grad.norm() > 1.e-12) + deallog << "Error grads" << std::endl; } } deallog << std::endl; @@ -77,6 +79,60 @@ check_support_points (const FiniteElement& fe) } +template +void +check_face_support_points (const FiniteElement& fe) +{ + deallog << "Face " << fe.get_name() << std::endl; + + const std::vector< Point > & sub_points = fe.get_unit_face_support_points(); + std::vector weights(sub_points.size()); + Quadrature sub_quadrature(sub_points, weights); + std::vector< Point > points(sub_points.size()); + + Triangulation tr; + GridGenerator::hyper_cube(tr); + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + MappingCartesian mapping; + + UpdateFlags flags = update_values | update_gradients; + FEFaceValues vals(mapping, fe, sub_quadrature, flags); + + for (unsigned int face=0;face::faces_per_cell;++face) + { + QProjector::project_to_face(sub_quadrature, face, points); + vals.reinit(dof.begin_active(), face); + + for (unsigned int k=0;k& p = points[k]; + deallog << p; + for (unsigned int i=0;i 1.e-12) + deallog << "Error values" << std::endl; + Tensor<1,dim> grad = fe.shape_grad_component(i,p,d); + grad -= vals.shape_grad_component(i,k,d); + if (grad.norm() > 1.e-12) + deallog << "Error grads" << std::endl; + } + } + deallog << std::endl; + } + deallog << std::endl; + } +} + + int main() { @@ -87,6 +143,8 @@ main() FE_RaviartThomasNodal<2> e20(0); check_support_points(e20); + check_face_support_points(e20); FE_RaviartThomasNodal<2> e21(1); check_support_points(e21); + check_face_support_points(e21); }