</p>
<ol>
+
+ <li> Changed: FEValues::normal_vector() for historical reasons returned a
+ value of type Point, though a normal vector is more adequately described
+ as a Tensor@<1,dim@>. Many similar cases were already clarified in deal.II
+ 8.3. The current case has now also been changed: FEValues::normal_vector()
+ now returns a Tensor, rather than a Point.
+ <br>
+ In a similar spirit, the FEValues::get_normal_vectors() function that
+ still returns a vector of Points has been deprecated and a new function,
+ FEValues::get_all_normal_vectors(), that returns a vector of tensors,
+ has been added. This was necessary since there is no way to change the
+ return type of the existing function in a backward compatible way. The
+ old function will be removed in the next version, and the new function
+ will then be renamed to the old name.
+ <br>
+ (Wolfgang Bangerth, 2015/08/20)
+ </li>
<li> Changed: The mesh_converter program has been removed from the
contrib folder. The equivalent functionality can now be found in
Vector<double> &local_vector = dinfo.vector(0).block(0);
const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();
std::vector<double> g(fe_v.n_quadrature_points);
// solution and the right hand side does not receive any contributions.
const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();
for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
{
/* ---------------------------------------------------------------------
*
- * Copyright (C) 2007 - 2014 by the deal.II authors
+ * Copyright (C) 2007 - 2015 by the deal.II authors
*
* This file is part of the deal.II library.
*
Vector<double> &cell_vector) const
{
const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
std::vector<double> g(fe_v.n_quadrature_points);
FullMatrix<double> &ue_ve_matrix) const
{
const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
// $\alpha$. It's form has also been given already in the introduction:
template <typename InputVector>
static
- void numerical_normal_flux (const Point<dim> &normal,
+ void numerical_normal_flux (const Tensor<1,dim> &normal,
const InputVector &Wplus,
const InputVector &Wminus,
const double alpha,
static
void
compute_Wminus (const BoundaryKind (&boundary_kind)[n_components],
- const Point<dim> &normal_vector,
+ const Tensor<1,dim> &normal_vector,
const DataVector &Wplus,
const Vector<double> &boundary_values,
const DataVector &Wminus)
fe_v.reinit(cell);
cell->get_dof_indices(local_dof_indices);
- const std::vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors();
+ const std::vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors();
wind.vector_value_list(q_points, cell_wind);
// We then form the integral over the current cell for all degrees of
std::vector<Vector<double> > singular_cell_wind( singular_quadrature.size(),
Vector<double>(dim) );
- const std::vector<Point<dim> > &singular_normals = fe_v_singular.get_normal_vectors();
- const std::vector<Point<dim> > &singular_q_points = fe_v_singular.get_quadrature_points();
+ const std::vector<Tensor<1,dim> > &singular_normals = fe_v_singular.get_all_normal_vectors();
+ const std::vector<Point<dim> > &singular_q_points = fe_v_singular.get_quadrature_points();
wind.vector_value_list(singular_q_points, singular_cell_wind);
{
fe_v.reinit(cell);
- const std::vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors();
+ const std::vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors();
cell->get_dof_indices(dofs);
fe_v.get_function_values(phi, local_phi);
const double JxW = scratch.fe_face_values.JxW(q);
const Point<dim> quadrature_point =
scratch.fe_face_values.quadrature_point(q);
- const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
+ const Tensor<1,dim> normal = scratch.fe_face_values.normal_vector(q);
const Tensor<1,dim> convection
= scratch.convection_velocity.value(quadrature_point);
std::vector<Point<spacedim> > quadrature_points;
/**
- * List of outward normal vectors at the quadrature points. This field is
- * filled in by the finite element class.
+ * List of outward normal vectors at the quadrature points.
*/
- std::vector<Point<spacedim> > normal_vectors;
+ std::vector<Tensor<1,spacedim> > normal_vectors;
/**
- * List of boundary forms at the quadrature points. This field is filled in
- * by the finite element class.
+ * List of boundary forms at the quadrature points.
*/
std::vector<Tensor<1,spacedim> > boundary_forms;
};
* For a face, return the outward normal vector to the cell at the
* <tt>i</tt>th quadrature point.
*
- * For a cell of codimension one, return the normal vector, as it is
- * specified by the numbering of the vertices.
+ * For a cell of codimension one, return the normal vector. There
+ * are of course two normal directions to a manifold in that case,
+ * and this function returns the "up" direction as induced by the
+ * numbering of the vertices.
*
* The length of the vector is normalized to one.
*
* @dealiiRequiresUpdateFlags{update_normal_vectors}
*/
- const Point<spacedim> &normal_vector (const unsigned int i) const;
+ const Tensor<1,spacedim> &normal_vector (const unsigned int i) const;
/**
* Return the normal vectors at the quadrature points. For a face, these are
* the orientation is given by the numbering of vertices.
*
* @dealiiRequiresUpdateFlags{update_normal_vectors}
+ *
+ * @note This function should really be named get_normal_vectors(),
+ * but this function already exists with a different return type
+ * that returns a vector of Point objects, rather than a vector of
+ * Tensor objects. This is a historical accident, but can not
+ * be fixed in a backward compatible style. That said, the
+ * get_normal_vectors() function is now deprecated, will be removed
+ * in the next version, and the current function will then be renamed.
*/
- const std::vector<Point<spacedim> > &get_normal_vectors () const;
+ const std::vector<Tensor<1,spacedim> > &get_all_normal_vectors () const;
+
+ /**
+ * Return the normal vectors at the quadrature points as a vector of
+ * Point objects. This function is deprecated because normal vectors
+ * are correctly represented by rank-1 Tensor objects, not Point objects.
+ * Use get_all_normal_vectors() instead.
+ *
+ * @dealiiRequiresUpdateFlags{update_normal_vectors}
+ *
+ * @deprecated
+ */
+ std::vector<Point<spacedim> > get_normal_vectors () const DEAL_II_DEPRECATED;
/**
* Transform a set of vectors, one for each quadrature point. The
template <int dim, int spacedim>
inline
-const Point<spacedim> &
+const Tensor<1,spacedim> &
FEValuesBase<dim,spacedim>::normal_vector (const unsigned int i) const
{
typedef FEValuesBase<dim,spacedim> FVB;
const CellSimilarity::Similarity cell_similarity,
const InternalData &data,
std::vector<Point<dim> > &quadrature_points,
- std::vector<Point<dim> > &normal_vectors) const;
+ std::vector<Tensor<1,dim> > &normal_vectors) const;
/**
* Value to indicate that a given face or subface number is invalid.
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = factor * fe.JxW(k);
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int d1=0; d1<dim; ++d1)
for (unsigned int d2=0; d2<dim; ++d2)
{
// v . nabla u n
- M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d1)[d2] *n(d2)* v;
+ M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d1)[d2] *n[d2]* v;
// v (nabla u)^T n
- M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d2)[d1] *n(d2)* v;
+ M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d2)[d1] *n[d2]* v;
// u nabla v n
- M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d1)[d2] *n(d2)* u;
+ M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d1)[d2] *n[d2]* u;
// u (nabla v)^T n
- M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d2)[d1] *n(d2)* u;
+ M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d2)[d1] *n[d2]* u;
}
}
}
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = factor * fe.JxW(k);
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int d1=0; d1<dim; ++d1)
{
for (unsigned int d2=0; d2<dim; ++d2)
{
// v . nabla u n
- result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
+ result(i) -= .5*dx* v * Dinput[d1][k][d2] * n[d2];
// v . (nabla u)^T n
- result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
+ result(i) -= .5*dx* v * Dinput[d2][k][d1] * n[d2];
// u nabla v n
- result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d1)[d2] * n(d2);
+ result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d1)[d2] * n[d2];
// u (nabla v)^T n
- result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d2)[d1] * n(d2);
+ result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d2)[d1] * n[d2];
}
}
}
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = factor * fe.JxW(k);
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int d1=0; d1<dim; ++d1)
{
for (unsigned int d2=0; d2<dim; ++d2)
{
// v . nabla u n
- result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
+ result(i) -= .5*dx* v * Dinput[d1][k][d2] * n[d2];
// v . (nabla u)^T n
- result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
+ result(i) -= .5*dx* v * Dinput[d2][k][d1] * n[d2];
// u nabla v n
- result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d1)[d2] * n(d2);
+ result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d1)[d2] * n[d2];
// u (nabla v)^T n
- result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d2)[d1] * n(d2);
+ result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d2)[d1] * n[d2];
}
}
}
for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
{
const double dx = fe1.JxW(k);
- const Point<dim> &n = fe1.normal_vector(k);
+ const Tensor<1,dim> n = fe1.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int d1=0; d1<dim; ++d1)
for (unsigned int d2=0; d2<dim; ++d2)
{
// v . nabla u n
- M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
- M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
- M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
- M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
+ M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n[d2] * v1;
+ M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n[d2] * v1;
+ M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n[d2] * v2;
+ M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n[d2] * v2;
// v (nabla u)^T n
- M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
- M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
- M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
- M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
+ M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n[d2] * v1;
+ M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n[d2] * v1;
+ M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n[d2] * v2;
+ M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n[d2] * v2;
// u nabla v n
- M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
- M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
- M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
- M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
+ M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n[d2] * u1;
+ M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n[d2] * u2;
+ M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n[d2] * u1;
+ M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n[d2] * u2;
// u (nabla v)^T n
- M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
- M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
- M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
- M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
+ M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n[d2] * u1;
+ M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n[d2] * u2;
+ M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n[d2] * u1;
+ M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n[d2] * u2;
}
}
}
for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
{
const double dx = fe1.JxW(k);
- const Point<dim> &n = fe1.normal_vector(k);
+ const Tensor<1,dim> n = fe1.normal_vector(k);
for (unsigned int i=0; i<n1; ++i)
for (unsigned int d1=0; d1<dim; ++d1)
for (unsigned int d2=0; d2<dim; ++d2)
{
// v . nabla u n
- result1(i) -= .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v1;
- result2(i) += .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v2;
+ result1(i) -= .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n[d2] * v1;
+ result2(i) += .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n[d2] * v2;
// v . (nabla u)^T n
- result1(i) -= .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v1;
- result2(i) += .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v2;
+ result1(i) -= .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n[d2] * v1;
+ result2(i) += .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n[d2] * v2;
// u nabla v n
- result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
- result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
+ result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d1)[d2] * n[d2] * (u1-u2);
+ result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d1)[d2] * n[d2] * (u1-u2);
// u (nabla v)^T n
- result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
- result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
+ result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d2)[d1] * n[d2] * (u1-u2);
+ result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d2)[d1] * n[d2] * (u1-u2);
}
}
}
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = fe.JxW(k) * factor;
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int d=0; d<n_comp; ++d)
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = factor * fe.JxW(k);
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
{
const double dnv = fe.shape_grad(i,k) * n;
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = factor * fe.JxW(k);
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int d=0; d<n_comp; ++d)
{
for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
{
const double dx = fe1.JxW(k);
- const Point<dim> &n = fe1.normal_vector(k);
+ const Tensor<1,dim> n = fe1.normal_vector(k);
for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
{
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
{
const double dx = fe1.JxW(k);
- const Point<dim> &n = fe1.normal_vector(k);
+ const Tensor<1,dim> n = fe1.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
{
for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int d=0; d<dim; ++d)
{
- u1dotn += n(d)*fe1.shape_value_component(j,k,d);
- v1dotn += n(d)*fe1.shape_value_component(i,k,d);
- u2dotn += n(d)*fe2.shape_value_component(j,k,d);
- v2dotn += n(d)*fe2.shape_value_component(i,k,d);
-
- ngradu1n += n*fe1.shape_grad_component(j,k,d)*n(d);
- ngradv1n += n*fe1.shape_grad_component(i,k,d)*n(d);
- ngradu2n += n*fe2.shape_grad_component(j,k,d)*n(d);
- ngradv2n += n*fe2.shape_grad_component(i,k,d)*n(d);
+ u1dotn += n[d]*fe1.shape_value_component(j,k,d);
+ v1dotn += n[d]*fe1.shape_value_component(i,k,d);
+ u2dotn += n[d]*fe2.shape_value_component(j,k,d);
+ v2dotn += n[d]*fe2.shape_value_component(i,k,d);
+
+ ngradu1n += n*fe1.shape_grad_component(j,k,d)*n[d];
+ ngradv1n += n*fe1.shape_grad_component(i,k,d)*n[d];
+ ngradu2n += n*fe2.shape_grad_component(j,k,d)*n[d];
+ ngradv2n += n*fe2.shape_grad_component(i,k,d)*n[d];
}
for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
{
- const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n(d);
- const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n(d);
+ const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n[d];
+ const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n[d];
- const double ve = fe2.shape_value_component(i,k,d)-v2dotn*n(d);
- const double dnve = n * fe2.shape_grad_component(i,k,d)-ngradv2n*n(d);
+ const double ve = fe2.shape_value_component(i,k,d)-v2dotn*n[d];
+ const double dnve = n * fe2.shape_grad_component(i,k,d)-ngradv2n*n[d];
- const double ui = fe1.shape_value_component(j,k,d)-u1dotn*n(d);
- const double dnui = n * fe1.shape_grad_component(j,k,d)-ngradu1n*n(d);
+ const double ui = fe1.shape_value_component(j,k,d)-u1dotn*n[d];
+ const double dnui = n * fe1.shape_grad_component(j,k,d)-ngradu1n*n[d];
- const double ue = fe2.shape_value_component(j,k,d)-u2dotn*n(d);
- const double dnue = n * fe2.shape_grad_component(j,k,d)-ngradu2n*n(d);
+ const double ue = fe2.shape_value_component(j,k,d)-u2dotn*n[d];
+ const double dnue = n * fe2.shape_grad_component(j,k,d)-ngradu2n*n[d];
M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
{
const double dx = fe1.JxW(k);
- const Point<dim> &n = fe1.normal_vector(k);
+ const Tensor<1,dim> n = fe1.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
{
for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
{
const double dx = fe1.JxW(k);
- const Point<dim> &n = fe1.normal_vector(k);
+ const Tensor<1,dim> n = fe1.normal_vector(k);
for (unsigned int i=0; i<n1; ++i)
for (unsigned int d=0; d<n_comp; ++d)
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = factor * fe.JxW(k);
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int d=0; d<d_max; ++d)
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
{
const double dx = factor * fe.JxW(k);
- const Point<dim> &n = fe.normal_vector(k);
+ const Tensor<1,dim> n = fe.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int d=0; d<d_max; ++d)
for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
{
const double dx = fe1.JxW(k);
- const Point<dim> &n = fe1.normal_vector(k);
+ const Tensor<1,dim> n = fe1.normal_vector(k);
for (unsigned int i=0; i<n_dofs; ++i)
for (unsigned int j=0; j<n_dofs; ++j)
for (unsigned int d=0; d<d_max; ++d)
get_quadrature_points() const;
/**
- * Return a vector of quadrature points in real space on the given cell.
+ * Return a vector of normal vectors in real space on the given cell.
* For compatibility with FEEvaluation, it returns tensors of vectorized
* arrays, even though all components are equal.
*/
// This avoids possible issues with the computation of the tangent.
// Store the normal at this quad point:
- Point<dim> normal_at_q_point = fe_face_values.normal_vector(q_point);
+ Tensor<1,dim> normal_at_q_point = fe_face_values.normal_vector(q_point);
for (unsigned int j = 0; j < associated_edge_dofs; ++j)
{
const unsigned int j_face_idx = associated_edge_dof_to_face_dof[j];
// Using n cross phi
edge_matrix(i,j)
+= fe_face_values.JxW (q_point)
- * ((phi_i[1]*normal_at_q_point(0) - phi_i[0]*normal_at_q_point(1))
- * (phi_j[1]*normal_at_q_point(0) - phi_j[0]*normal_at_q_point(1)));
+ * ((phi_i[1]*normal_at_q_point[0] - phi_i[0]*normal_at_q_point[1])
+ * (phi_j[1]*normal_at_q_point[0] - phi_j[0]*normal_at_q_point[1]));
}
// Using n cross phi
edge_rhs(j)
+= fe_face_values.JxW (q_point)
- * ((values[q_point] (first_vector_component+1) * normal_at_q_point (0)
- - values[q_point] (first_vector_component) * normal_at_q_point (1))
- * (phi_j[1]*normal_at_q_point(0) - phi_j[0]*normal_at_q_point(1)));
+ * ((values[q_point] (first_vector_component+1) * normal_at_q_point[0]
+ - values[q_point] (first_vector_component) * normal_at_q_point[1])
+ * (phi_j[1]*normal_at_q_point[0] - phi_j[0]*normal_at_q_point[1]));
}
}
// associated with this face. We also must include the residuals from the
// shape funcations associated with edges.
Tensor<1, dim> tmp;
- Tensor<1, dim> normal_vector,
- cross_product_i,
+ Tensor<1, dim> cross_product_i,
cross_product_j,
cross_product_rhs;
- // Store all normal vectors at quad points:
- std::vector<Point<dim> > normal_vector_list(fe_face_values.get_normal_vectors());
-
// Loop to construct face linear system.
for (unsigned int q_point = 0;
q_point < fe_face_values.n_quadrature_points; ++q_point)
}
// Tensor of normal vector on the face at q_point;
- for (unsigned int d = 0; d < dim; ++d)
- {
- normal_vector[d] = normal_vector_list[q_point](d);
- }
+ const Tensor<1,dim> normal_vector = fe_face_values.normal_vector(q_point);
+
// Now compute the linear system:
// On a face:
// The matrix entries are:
// functions supported on the boundary.
const FEValuesExtractors::Vector vec (first_vector_component);
const FiniteElement<2> &fe = cell->get_fe ();
- const std::vector<Point<2> > &normals = fe_values.get_normal_vectors ();
+ const std::vector<Tensor<1,2> > &normals = fe_values.get_all_normal_vectors ();
const unsigned int
face_coordinate_direction[GeometryInfo<2>::faces_per_cell] = {1, 1, 0, 0};
std::vector<Vector<double> >
// functions supported on the boundary.
const FEValuesExtractors::Vector vec (first_vector_component);
const FiniteElement<3> &fe = cell->get_fe ();
- const std::vector<Point<3> > &normals = fe_values.get_normal_vectors ();
+ const std::vector<Tensor<1,3> > &normals = fe_values.get_all_normal_vectors ();
const unsigned int
face_coordinate_directions[GeometryInfo<3>::faces_per_cell][2] = {{1, 2},
{1, 2},
= (cell->face(face_no)->get_boundary().normal_vector
(cell->face(face_no),
fe_values.quadrature_point(i)));
- if (normal_vector * static_cast<Tensor<1,dim> >(fe_values.normal_vector(i)) < 0)
+ if (normal_vector * fe_values.normal_vector(i) < 0)
normal_vector *= -1;
Assert (std::fabs(normal_vector.norm() - 1) < 1e-14,
ExcInternalError());
{
// compute (f.n) n
const Number f_dot_n
- = (data.psi_grads[q][k] * Tensor<1,spacedim,Number>(fe_values.normal_vector(q)));
- const Tensor<1,spacedim,Number> f_dot_n_times_n
- = f_dot_n * Tensor<1,spacedim,Number>(fe_values.normal_vector(q));
+ = (data.psi_grads[q][k] * fe_values.normal_vector(q));
+ const Tensor<1,spacedim,Number> f_dot_n_times_n (f_dot_n * fe_values.normal_vector(q));
data.psi_grads[q][k] -= (data.function_grads[q][k] + f_dot_n_times_n);
}
template <int dim, int spacedim>
-const std::vector<Point<spacedim> > &
-FEValuesBase<dim,spacedim>::get_normal_vectors () const
+const std::vector<Tensor<1,spacedim> > &
+FEValuesBase<dim,spacedim>::get_all_normal_vectors () const
{
typedef FEValuesBase<dim,spacedim> FEVB;
Assert (this->update_flags & update_normal_vectors,
+template <int dim, int spacedim>
+std::vector<Point<spacedim> >
+FEValuesBase<dim,spacedim>::get_normal_vectors () const
+{
+ typedef FEValuesBase<dim,spacedim> FEVB;
+ Assert (this->update_flags & update_normal_vectors,
+ typename FEVB::ExcAccessToUninitializedField("update_normal_vectors"));
+
+ // copy things into a vector of Points, then return that
+ std::vector<Point<spacedim> > tmp (this->normal_vectors.size());
+ for (unsigned int q=0; q<this->normal_vectors.size(); ++q)
+ tmp[q] = Point<spacedim>(this->normal_vectors[q]);
+
+ return tmp;
+}
+
+
+
template <int dim, int spacedim>
void
FEValuesBase<dim,spacedim>::transform (
const CellSimilarity::Similarity cell_similarity,
const InternalData &data,
std::vector<Point<dim> > &quadrature_points,
- std::vector<Point<dim> > &normal_vectors) const
+ std::vector<Tensor<1,dim> > &normal_vectors) const
{
const UpdateFlags update_flags = data.update_each;
Assert (dynamic_cast<const InternalData *> (&internal_data) != 0, ExcInternalError());
const InternalData &data = static_cast<const InternalData &> (internal_data);
- std::vector<Point<dim> > dummy;
+ std::vector<Tensor<1,dim> > dummy;
compute_fill (cell, invalid_face_number, invalid_face_number, cell_similarity,
data,
// thus does not depend on the number of components of the data
// vector
if (update_flags & update_normal_vectors)
- data.patch_normals=this_fe_patch_values.get_normal_vectors();
+ {
+//TODO: undo this copying when we can change the data type of
+// data.patch_normals to Tensor<1,spacedim> as well
+ for (unsigned int q=0; q<this_fe_patch_values.n_quadrature_points; ++q)
+ data.patch_normals[q] = Point<dim>(this_fe_patch_values.get_all_normal_vectors()[q]);
+ }
if (n_components == 1)
{
/**
* The normal vectors of the finite element function on one face
*/
- std::vector<Point<spacedim> > normal_vectors;
+ std::vector<Tensor<1,spacedim> > normal_vectors;
/**
* Two arrays needed for the values of coefficients in the jumps, if
// change the sign. We take the outward normal.
parallel_data.normal_vectors =
- fe_face_values_cell.get_present_fe_values().get_normal_vectors();
+ fe_face_values_cell.get_present_fe_values().get_all_normal_vectors();
for (unsigned int n=0; n<n_solution_vectors; ++n)
for (unsigned int component=0; component<n_components; ++component)
if (!base.conforms(FiniteElementData<dim>::H1) &&
base.conforms(FiniteElementData<dim>::Hdiv))
for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- normal_adjustment[point][comp] = fe_values.normal_vector(point)(bcomp)
- * fe_values.normal_vector(point)(bcomp);
+ normal_adjustment[point][comp] = fe_values.normal_vector(point)[bcomp]
+ * fe_values.normal_vector(point)[bcomp];
}
for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
fe_face_values2.JxW(q)) < 1e-14,
ExcInternalError());
AssertThrow ((fe_face_values1.normal_vector(q) +
- fe_face_values2.normal_vector(q)).square()
+ fe_face_values2.normal_vector(q)).norm_square()
< 1e-20,
ExcInternalError());
}
for (unsigned int q_point=0; q_point<n_q_face; ++q_point)
{
- Point<2> vn = fe_v_face.normal_vector (q_point);
- double nx = vn(0);
- double ny = vn(1);
+ Tensor<1,2> vn = fe_v_face.normal_vector (q_point);
+ double nx = vn[0];
+ double ny = vn[1];
double u = this_value[q_point](0);
double v = this_value[q_point](1);
// ---------------------------------------------------------------------
//
-// Copyright (C) 2005 - 2014 by the deal.II authors
+// Copyright (C) 2005 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
Vector<double> &cell_vector) const
{
const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
std::vector<double> g(fe_v.n_quadrature_points);
FullMatrix<double> &ue_vi_matrix) const
{
const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
beta_function.value_list (fe_v.get_quadrature_points(), beta);
FullMatrix<double> &ue_ve_matrix) const
{
const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
const double JxW = scratch.fe_face_values.JxW(q);
const Point<dim> quadrature_point =
scratch.fe_face_values.quadrature_point(q);
- const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
+ const Tensor<1,dim> normal = scratch.fe_face_values.normal_vector(q);
const Tensor<1,dim> convection
= scratch.convection_velocity.value(quadrature_point);
const double JxW = scratch.fe_face_values.JxW(q);
const Point<dim> quadrature_point =
scratch.fe_face_values.quadrature_point(q);
- const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
+ const Tensor<1,dim> normal = scratch.fe_face_values.normal_vector(q);
const Tensor<1,dim> convection
= scratch.convection_velocity.value(quadrature_point);
FullMatrix<double> cell_DLP_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
- std::vector< Point<spacedim> > cell_normals_i, cell_normals_j;
+ std::vector< Tensor<1,spacedim> > cell_normals_i, cell_normals_j;
std::vector<types::global_dof_index> local_dof_indices_i (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices_j (dofs_per_cell);
{
fe_values_i.reinit (cell_i);
- cell_normals_i = fe_values_i.get_normal_vectors();
+ cell_normals_i = fe_values_i.get_all_normal_vectors();
cell_i->get_dof_indices (local_dof_indices_i);
// if (cell_i->index()%100==0)
for (cell_j=dof_handler.begin_active(); cell_j!=endc; ++cell_j)
{
fe_values_j.reinit (cell_j);
- cell_normals_j = fe_values_j.get_normal_vectors();
+ cell_normals_j = fe_values_j.get_all_normal_vectors();
cell_j->get_dof_indices (local_dof_indices_j);
if (cell_j != cell_i)
update_gradients |
update_normal_vectors );
- std::vector< Point<spacedim> > cell_normals(q_iterated.size());
+ std::vector< Tensor<1,spacedim> > cell_normals(q_iterated.size());
std::vector< Point<spacedim> > cell_tangentials(q_iterated.size());
std::vector<types::global_dof_index> local_dof_indices (fe_q.dofs_per_cell);
// std::cout<<std::endl;
- cell_normals = fe_values_q.get_normal_vectors();
+ cell_normals = fe_values_q.get_all_normal_vectors();
for (unsigned int i=0; i<q_iterated.size(); ++i)
{
cell_tangentials[i][0] = cell_normals[i][1];
double term_S(const Tensor<1,3> &r,
const Tensor<1,3> &a1,
const Tensor<1,3> &a2,
- const Point<3> &n,
+ const Tensor<1,3> &n,
const double &rn_c);
double term_D(const Tensor<1,3> &r,
ExcDimensionMismatch(dst.size(), 2));
fe_values->reinit(cell);
vector<DerivativeForm<1,2,3> > jacobians = fe_values->get_jacobians();
- vector<Point<3> > normals = fe_values->get_normal_vectors();
+ vector<Tensor<1,3> > normals = fe_values->get_all_normal_vectors();
- Point<3> n,n_c;
+ Tensor<1,3> n,n_c;
Tensor<1,3> r_c = point-cell->center();
n_c = normals[4];
LaplaceKernelIntegration<dim>::term_S (const Tensor<1,3> &r,
const Tensor<1,3> &a1,
const Tensor<1,3> &a2,
- const Point<3> &n,
+ const Tensor<1,3> &n,
const double &rn_c)
{
Point<3> ra1, ra2, a12;
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- std::vector< Point<spacedim> > cell_normals(q_midpoint.size());
+ std::vector< Tensor<1,spacedim> > cell_normals(q_midpoint.size());
std::vector< Point<spacedim> > cell_tangentials(q_midpoint.size());
std::vector<double> shape_directional_derivative(dofs_per_cell);
Vector<double> projected_directional_derivative(triangulation.n_cells());
fe_values.reinit(cell);
cell-> get_dof_indices (local_dof_indices);
- cell_normals = fe_values.get_normal_vectors();
+ cell_normals = fe_values.get_all_normal_vectors();
// The cell tangential is calculated
// in the midpoint of the cell. For
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
- const std::vector<Point<spacedim> > &cellnormals = fe_values.get_normal_vectors();
+ const std::vector<Tensor<1,spacedim> > &cellnormals = fe_values.get_all_normal_vectors();
const std::vector<Point<spacedim> > &quad_points = fe_values.get_quadrature_points();
for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
// Neumann storage
std::vector<Vector<double> > neumann_value_list(n_face_q_points, Vector<double>(fe.n_components()));
- std::vector<Point<dim> > normal_vector_list(fe_face_values.get_normal_vectors());
+ std::vector<Tensor<1,dim> > normal_vector_list(fe_face_values.get_all_normal_vectors());
Tensor<1,dim> neumann_value_vector(dim);
Tensor<1,dim> neumann_value(dim);
Tensor<1,dim> normal_vector;
// ---------------------------------------------------------------------
//
-// Copyright (C) 2007 - 2014 by the deal.II authors
+// Copyright (C) 2007 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
{
double diff = 0.;
for (unsigned int d=0; d<dim; ++d)
- diff += fe.normal_vector(k)(d) * (f_values[k](d) - fe_values[k](d));
+ diff += fe.normal_vector(k)[d] * (f_values[k](d) - fe_values[k](d));
result += fe.JxW(k) * diff * diff;
}
}
// ---------------------------------------------------------------------
//
-// Copyright (C) 1998 - 2014 by the deal.II authors
+// Copyright (C) 1998 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
JumpFunction () : Function<dim>(1) {}
double jump (const Point< dim> &p,
- const Point<dim> &normal) const;
+ const Tensor<1,dim> &normal) const;
};
template <int dim>
double
JumpFunction<dim>::jump (const Point<dim> &p,
- const Point<dim> &normal) const
+ const Tensor<1,dim> &normal) const
{
double x = p[0];
double y = p[1];
{
fe_v_face.reinit (cell, face_no);
- const std::vector<Point<dim> > &normals =
- fe_v_face.get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals =
+ fe_v_face.get_all_normal_vectors ();
const std::vector<Point<dim> > &quad_points=
fe_v_face.get_quadrature_points();
for (unsigned int nx=0; nx<nq; ++nx)
{
const Point<dim> x = fe_values.quadrature_point(k);
- const Point<dim> &n = fe_values.normal_vector(k);
+ const Tensor<1,dim> n = fe_values.normal_vector(k);
const double ds = fe_values.JxW(k);
deallog << x << '\t' << n << '\t' << ds << std::endl;
{
fe_values.reinit(cell, face_nr, sub_nr);
- const std::vector<Point<dim> > &normals
- =fe_values.get_normal_vectors();
+ const std::vector<Tensor<1,dim> > &normals
+ =fe_values.get_all_normal_vectors();
unsigned int k=0;
for (unsigned int ny=0; ny<((dim>2) ? nq : 1); ++ny)
// there should now be two
// normal vectors, one for
// each vertex of the face
- Assert (c1_values.get_normal_vectors().size() == 2,
+ Assert (c1_values.get_all_normal_vectors().size() == 2,
ExcInternalError());
// check that these two
for (unsigned int q_point=0; q_point<n_q_face; ++q_point)
{
- Point<2> vn = fe_v_face.normal_vector (q_point);
- double nx = vn(0);
- double ny = vn(1);
+ Tensor<1,2> vn = fe_v_face.normal_vector (q_point);
+ double nx = vn[0];
+ double ny = vn[1];
double u = this_value[q_point + offset](0);
double v = this_value[q_point + offset](1);
for (unsigned int q_point=0; q_point<n_q_face; ++q_point)
{
- Point<2> vn = fe_v_face.normal_vector (q_point);
- double nx = vn(0);
- double ny = vn(1);
+ Tensor<1,2> vn = fe_v_face.normal_vector (q_point);
+ double nx = vn[0];
+ double ny = vn[1];
double u = this_value[q_point](0);
double v = this_value[q_point](1);
// ---------------------------------------------------------------------
//
-// Copyright (C) 2003 - 2014 by the deal.II authors
+// Copyright (C) 2003 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
// so their sum should be
// close to zero
Assert ((fe_face_values1.normal_vector(0) +
- fe_face_values2.normal_vector(0)).square()
+ fe_face_values2.normal_vector(0)).norm_square()
< 1e-20,
ExcInternalError());
}
// ---------------------------------------------------------------------
//
-// Copyright (C) 2005 - 2014 by the deal.II authors
+// Copyright (C) 2005 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
Vector<double> &cell_vector) const
{
const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_present_fe_values().get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
std::vector<double> g(fe_v.get_present_fe_values().n_quadrature_points);
FullMatrix<double> &ue_vi_matrix) const
{
const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_present_fe_values().get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
FullMatrix<double> &ue_ve_matrix) const
{
const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
+ const std::vector<Tensor<1,dim> > &normals = fe_v.get_present_fe_values().get_all_normal_vectors ();
std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
projects / dealii.git / commitdiff