// reference_function(x_j)-\psi_j
// and assign that to the vector
// \psi.
- const unsigned int n_dofs = fe.total_dofs;
const unsigned int n_q_points = q.n_quadrature_points;
- const dFMatrix & shape_values = fe_values.get_shape_values();
- vector<double> dof_values(fe.total_dofs, 0);
- cell->get_dof_values (solution, dof_values);
-
vector<double> psi;
// in praxi: first compute
psi);
// then subtract finite element
// solution
- for (unsigned int j=0; j<n_q_points; ++j)
- for (unsigned int i=0; i<n_dofs; ++i)
- psi[j] -= dof_values[i]*shape_values(i,j);
+ if (true)
+ {
+ vector<double> function_values (n_q_points, 0);
+ fe_values.get_function_values (solution, function_values);
+
+ transform (psi.begin(), psi.end(),
+ function_values.begin(),
+ psi.begin(),
+ minus<double>());
+ };
+
// for L1_norm and Linfty_norm:
// take absolute
// same procedure as above, but now
// psi is a vector of gradients
- const unsigned int n_dofs = fe.total_dofs;
const unsigned int n_q_points = q.n_quadrature_points;
- const vector<vector<Point<dim> > > & shape_grads = fe_values.get_shape_grads();
- vector<double> dof_values(fe.total_dofs, 0);
- cell->get_dof_values (solution, dof_values);
-
vector<Point<dim> > psi;
// in praxi: first compute
// then subtract finite element
// solution
- for (unsigned int j=0; j<n_q_points; ++j)
- for (unsigned int i=0; i<n_dofs; ++i)
- psi[j] -= dof_values[i]*shape_grads[i][j];
+ if (true)
+ {
+ vector<Point<dim> > function_grads (n_q_points, Point<dim>());
+ fe_values.get_function_grads (solution, function_grads);
+ transform (psi.begin(), psi.end(),
+ function_grads.begin(),
+ psi.begin(),
+ minus<Point<dim> >());
+ };
// take square of integrand
vector<double> psi_square (psi.size(), 0.0);
for (unsigned int i=0; i<n_q_points; ++i)
// number of integration points per face
const unsigned int n_q_points = quadrature.n_quadrature_points;
- // number of dofs per cell
- const unsigned int n_dofs = fe.total_dofs;
// make up a fe face values object for the
// restriction of the finite element function
// let psi be a short name for
// [grad u_h]
vector<Point<dim> > psi(n_q_points);
+ fe_face_values_cell.get_function_grads (solution, psi);
- // get a list of the values of
- // the solution for the ansatz
- // functions on this cell
- vector<double> dof_values(fe.total_dofs, 0);
- cell->get_dof_values (solution, dof_values);
-
- // get a list of the gradients of
- // the ansatz functions on this
- // cell at the quadrature points
- const vector<vector<Point<dim> > > &shape_grads(fe_face_values_cell.
- get_shape_grads());
-
- // compute the gradients of the solution
- // at the quadrature points by summing
- // over the ansatz functions.
- for (unsigned int j=0; j<n_q_points; ++j)
- for (unsigned int i=0; i<n_dofs; ++i)
- psi[j] += dof_values[i]*shape_grads[i][j];
-
-
// now compute over the other side of
// the face
fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor,
fe, boundary);
- // get a list of the values of
- // the solution for the ansatz
- // functions on this neighbor
- neighbor->get_dof_values (solution, dof_values);
- // get a list of the gradients of the
- //
- const vector<vector<Point<dim> > > &neighbor_grads (fe_face_values_cell.
- get_shape_grads());
- // subtract the gradients of the
- // solution on the neigbor cell
- // at the quadrature points from
- // those of the present cell
- for (unsigned int j=0; j<n_q_points; ++j)
- for (unsigned int i=0; i<n_dofs; ++i)
- psi[j] -= dof_values[i]*neighbor_grads[i][j];
+ // get a list of the gradients of
+ // the finite element solution
+ // restricted to the neighbor cell
+ vector<Point<dim> > neighbor_psi (n_q_points);
+ fe_face_values_neighbor.get_function_grads (solution, neighbor_psi);
+
+ // compute the jump in the gradients
+ transform (psi.begin(), psi.end(),
+ neighbor_psi.begin(),
+ psi.begin(),
+ minus<Point<dim> >());
};
vector<double> phi(n_q_points,0);
const vector<Point<dim> > &normal_vectors(fe_face_values_cell.
get_normal_vectors());
+
for (unsigned int point=0; point<n_q_points; ++point)
phi[point] = psi[point]*normal_vectors[point];