if (inverse_node_matrix.n_cols() == 0)
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
data->shape_values[i][k] = values[i];
- else
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- for (unsigned int j=0; j<this->dofs_per_cell; ++j)
- data->shape_values[i][k] = inverse_node_matrix(j,i) * values[j];
+ else
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<this->dofs_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)
else
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
for (unsigned int j=0; j<this->dofs_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;
}
for (unsigned int i=0; i<n_dofs; ++i)
{
- const unsigned int first = i*dim;//data.shape_function_to_row_table[i];
+ const unsigned int first = data.shape_function_to_row_table[i];
if (flags & update_values)
for (unsigned int k=0; k<n_quad; ++k)
MappingCartesian<dim> mapping;
Quadrature<dim> support_quadrature(points, weights);
- UpdateFlags flags = update_values;
+ UpdateFlags flags = update_values | update_gradients;
FEValues<dim> vals(mapping, fe, support_quadrature, flags);
vals.reinit(dof.begin_active());
for (unsigned int d=0;d<dim;++d)
{
const double sf = fe.shape_value_component(i,p,d);
- const double sv = vals.shape_value_component(i,k,d);
- deallog << ' ' << (int) round(60*sf)
- << '/' << (int) round(60*sv);
+ deallog << ' ' << (int) round(60*sf);
-// const double diff = std::abs(sf - vals.shape_value_component(i,k,d));
-// if (diff > 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;
}
+template <int dim>
+void
+check_face_support_points (const FiniteElement<dim>& fe)
+{
+ deallog << "Face " << fe.get_name() << std::endl;
+
+ const std::vector< Point<dim-1> > & sub_points = fe.get_unit_face_support_points();
+ std::vector<double> weights(sub_points.size());
+ Quadrature<dim-1> sub_quadrature(sub_points, weights);
+ std::vector< Point<dim> > points(sub_points.size());
+
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr);
+ DoFHandler<dim> dof(tr);
+ dof.distribute_dofs(fe);
+
+ MappingCartesian<dim> mapping;
+
+ UpdateFlags flags = update_values | update_gradients;
+ FEFaceValues<dim> vals(mapping, fe, sub_quadrature, flags);
+
+ for (unsigned int face=0;face<GeometryInfo<dim>::faces_per_cell;++face)
+ {
+ QProjector<dim>::project_to_face(sub_quadrature, face, points);
+ vals.reinit(dof.begin_active(), face);
+
+ for (unsigned int k=0;k<points.size();++k)
+ {
+ const Point<dim>& p = points[k];
+ deallog << p;
+ for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+ {
+ deallog << '\t';
+ for (unsigned int d=0;d<dim;++d)
+ {
+ const double sf = fe.shape_value_component(i,p,d);
+ deallog << ' ' << (int) round(60*sf);
+
+ 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;
+ }
+ deallog << std::endl;
+ }
+}
+
+
int
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);
}
--- /dev/null
+
+DEAL::FE_RaviartThomasNodal<2>(0)
+DEAL::0.500000 0.00000 0 -60 30 0 0 0 -30 0
+DEAL::1.00000 0.500000 0 -30 60 0 0 30 0 0
+DEAL::0.500000 1.00000 0 0 30 0 0 60 -30 0
+DEAL::0.00000 0.500000 0 -30 0 0 0 30 -60 0
+DEAL::
+DEAL::Face FE_RaviartThomasNodal<2>(0)
+DEAL::0.500000 0.00000 0 -60 30 0 0 0 -30 0
+DEAL::
+DEAL::1.00000 0.500000 0 -30 60 0 0 30 0 0
+DEAL::
+DEAL::0.500000 1.00000 0 0 30 0 0 60 -30 0
+DEAL::
+DEAL::0.00000 0.500000 0 -30 0 0 0 30 -60 0
+DEAL::
+DEAL::FE_RaviartThomasNodal<2>(1)
+DEAL::0.211325 0.00000 0 -60 0 0 -10 0 3 0 0 0 0 0 -37 0 10 0 55 0 -15 0 0 0 0 0
+DEAL::0.788675 0.00000 0 0 0 -60 37 0 -10 0 0 0 0 0 10 0 -3 0 55 0 -15 0 0 0 0 0
+DEAL::1.00000 0.211325 0 10 0 -37 60 0 0 0 0 3 0 -10 0 0 0 0 0 0 0 0 0 -15 0 55
+DEAL::1.00000 0.788675 0 -3 0 10 0 0 60 0 0 -10 0 37 0 0 0 0 0 0 0 0 0 -15 0 55
+DEAL::0.211325 1.00000 0 0 0 0 3 0 -10 0 0 60 0 0 10 0 -37 0 -15 0 55 0 0 0 0 0
+DEAL::0.788675 1.00000 0 0 0 0 -10 0 37 0 0 0 0 60 -3 0 10 0 -15 0 55 0 0 0 0 0
+DEAL::0.00000 0.211325 0 -37 0 10 0 0 0 0 0 -10 0 3 -60 0 0 0 0 0 0 0 0 55 0 -15
+DEAL::0.00000 0.788675 0 10 0 -3 0 0 0 0 0 37 0 -10 0 0 -60 0 0 0 0 0 0 55 0 -15
+DEAL::0.500000 0.211325 0 -14 0 -14 0 0 0 0 0 -4 0 -4 0 0 0 0 60 0 0 0 0 20 0 20
+DEAL::0.500000 0.788675 0 4 0 4 0 0 0 0 0 14 0 14 0 0 0 0 0 0 60 0 0 20 0 20
+DEAL::0.211325 0.500000 0 0 0 0 -4 0 -4 0 0 0 0 0 -14 0 -14 0 20 0 20 0 0 60 0 0
+DEAL::0.788675 0.500000 0 0 0 0 14 0 14 0 0 0 0 0 4 0 4 0 20 0 20 0 0 0 0 60
+DEAL::
+DEAL::Face FE_RaviartThomasNodal<2>(1)
+DEAL::0.211325 0.00000 0 -60 0 0 -10 0 3 0 0 0 0 0 -37 0 10 0 55 0 -15 0 0 0 0 0
+DEAL::0.788675 0.00000 0 0 0 -60 37 0 -10 0 0 0 0 0 10 0 -3 0 55 0 -15 0 0 0 0 0
+DEAL::
+DEAL::1.00000 0.211325 0 10 0 -37 60 0 0 0 0 3 0 -10 0 0 0 0 0 0 0 0 0 -15 0 55
+DEAL::1.00000 0.788675 0 -3 0 10 0 0 60 0 0 -10 0 37 0 0 0 0 0 0 0 0 0 -15 0 55
+DEAL::
+DEAL::0.211325 1.00000 0 0 0 0 3 0 -10 0 0 60 0 0 10 0 -37 0 -15 0 55 0 0 0 0 0
+DEAL::0.788675 1.00000 0 0 0 0 -10 0 37 0 0 0 0 60 -3 0 10 0 -15 0 55 0 0 0 0 0
+DEAL::
+DEAL::0.00000 0.211325 0 -37 0 10 0 0 0 0 0 -10 0 3 -60 0 0 0 0 0 0 0 0 55 0 -15
+DEAL::0.00000 0.788675 0 10 0 -3 0 0 0 0 0 37 0 -10 0 0 -60 0 0 0 0 0 0 55 0 -15
+DEAL::