* #i#th shape function for the mapping
* from unit to real cell.
*/
- virtual Point<dim> shape_grad_transform (const unsigned int i,
- const Point<dim> &p) const = 0;
+ virtual Tensor<1,dim> shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const = 0;
/**
* Compute the Jacobian matrix and the
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &shape_grads_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
const Boundary<dim> &boundary) const;
/**
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &shape_grads_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
const Boundary<dim> &boundary) const;
/**
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &shape_grads_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_grads_transform,
const Boundary<dim> &boundary) const;
/**
* information on this function.
*/
virtual void get_support_points (const DoFHandler<dim>::cell_iterator &cell,
- const Boundary<dim> &boundary,
- vector<Point<dim> > &support_points) const;
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* information on this function.
*/
virtual void get_face_support_points (const DoFHandler<dim>::face_iterator &face,
- const Boundary<dim> &boundary,
- vector<Point<dim> > &support_points) const;
+ const Boundary<dim> &boundary,
+ vector<Point<dim> > &support_points) const;
/**
* Refer to the base class for detailed
* computation of the transformation from
* unit cell to real space cell.
*/
- virtual Point<dim> shape_grad_transform (const unsigned int i,
- const Point<dim> &p) const;
+ virtual Tensor<1,dim> shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &shape_grad_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
const Boundary<dim> &boundary) const;
DeclException0 (ExcNotUseful);
* computation of the transformation from
* unit cell to real space cell.
*/
- virtual Point<dim> shape_grad_transform (const unsigned int i,
- const Point<dim> &p) const;
+ virtual Tensor<1,dim> shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &shape_grad_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
const Boundary<dim> &boundary) const;
};
* function at the #j# quadrature point.
* If you want to get the derivative in
* one of the coordinate directions, use
- * the appropriate function of the #Point#
+ * the appropriate function of the #Tensor#
* class to extract one component. Since
* only a reference to the gradient's value
* is returned, there should be no major
* #gradients# object already has the
* right size.
*/
- void get_function_grads (const dVector &fe_function,
- vector<Point<dim> > &gradients) const;
+ void get_function_grads (const dVector &fe_function,
+ vector<Tensor<1,dim> > &gradients) const;
/**
* Return the position of the #i#th
* Analogous to the #shape_values_transform#
* array of the base class.
*/
- vector<vector<Point<dim> > > unit_shape_gradients_transform;
+ vector<vector<Tensor<1,dim> > > unit_shape_gradients_transform;
/**
* Array of quadrature points in the unit
* Analogous to the #shape_values_transform#
* array of the base class.
*/
- vector<vector<vector<Point<dim> > > > unit_shape_gradients_transform;
+ vector<vector<vector<Tensor<1,dim> > > > unit_shape_gradients_transform;
/**
* Array of quadrature points on the
#include <vector>
#include <base/point.h>
#include <grid/geometry_info.h>
-
+#include <base/subscriptor.h>
//forward declaration needed
* @author Wolfgang Bangerth, 1998
*/
template <int dim>
-class Triangulation : public TriaDimensionInfo<dim> {
+class Triangulation : public TriaDimensionInfo<dim>, public Subscriptor {
public:
// insert these definitions for gcc2.8,
// since it can't inherit typedefs (I
template <int dim>
DoFHandler<dim>::DoFHandler (Triangulation<dim> *tria) :
tria(tria),
- used_dofs (0) {};
+ used_dofs (0)
+{
+ tria->subscribe ();
+};
template <int dim>
-DoFHandler<dim>::~DoFHandler ()
-{}
+DoFHandler<dim>::~DoFHandler () {
+ tria->unsubscribe ();
+}
vector<Point<1> > &q_points,
const bool compute_q_points,
const dFMatrix &,
- const vector<vector<Point<1> > > &,
+ const vector<vector<Tensor<1,1> > > &,
const Boundary<1> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
vector<Point<1> > &,
const bool,
const dFMatrix &,
- const vector<vector<Point<1> > > &,
+ const vector<vector<Tensor<1,1> > > &,
const Boundary<1> &) const {
Assert (false, ExcNotImplemented());
};
vector<Point<1> > &,
const bool,
const dFMatrix &,
- const vector<vector<Point<1> > > &,
+ const vector<vector<Tensor<1,1> > > &,
const Boundary<1> &) const {
Assert (false, ExcNotImplemented());
};
vector<Point<dim> > &,
const bool,
const dFMatrix &,
- const vector<vector<Point<dim> > > &,
+ const vector<vector<Tensor<1,dim> > > &,
const Boundary<dim> &) const {
Assert (false, ExcPureFunctionCalled());
};
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &shape_gradients_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_gradients_transform,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &shape_gradients_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_gradients_transform,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
template <>
void FECrissCross<1>::get_support_points (const DoFHandler<1>::cell_iterator &,
- const Boundary<1> &,
- vector<Point<1> > &) const {
+ const Boundary<1> &,
+ vector<Point<1> > &) const {
Assert (false, ExcNotUseful());
};
template <>
void FECrissCross<1>::get_face_support_points (const DoFHandler<1>::face_iterator &,
- const Boundary<1> &,
- vector<Point<1> > &) const {
+ const Boundary<1> &,
+ vector<Point<1> > &) const {
Assert (false, ExcNotUseful());
};
template <>
-Point<1> FECrissCross<1>::shape_grad_transform (const unsigned int,
- const Point<1> &) const {
+Tensor<1,1> FECrissCross<1>::shape_grad_transform (const unsigned int,
+ const Point<1> &) const {
Assert (false, ExcNotUseful());
return Point<1>();
};
vector<Point<1> > &,
const bool,
const dFMatrix &,
- const vector<vector<Point<1> > > &,
+ const vector<vector<Tensor<1,1> > > &,
const Boundary<1> &) const {
Assert (false, ExcNotUseful());
};
template <>
void FECrissCross<2>::get_support_points (const DoFHandler<2>::cell_iterator &cell,
- const Boundary<2> &,
- vector<Point<2> > &support_points) const {
+ const Boundary<2> &,
+ vector<Point<2> > &support_points) const {
const unsigned int dim = 2;
Assert (support_points.size() == total_dofs,
template <>
void FECrissCross<2>::get_face_support_points (const DoFHandler<2>::face_iterator &face,
- const Boundary<2> &,
- vector<Point<2> > &support_points) const {
+ const Boundary<2> &,
+ vector<Point<2> > &support_points) const {
const unsigned int dim = 2;
Assert ((support_points.size() == dofs_per_face) &&
template <>
inline
double FECrissCross<2>::shape_value_transform (const unsigned int i,
- const Point<2> &p) const {
+ const Point<2> &p) const {
// use an isoparametric ansatz
return shape_value(i,p);
};
template <>
-Point<2> FECrissCross<2>::shape_grad_transform (const unsigned int i,
- const Point<2> &p) const {
+Tensor<1,2> FECrissCross<2>::shape_grad_transform (const unsigned int i,
+ const Point<2> &p) const {
// use an isoparametric ansatz
return shape_grad(i,p);
};
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &/*shape_grad_transform*/,
+ const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
vertices[l] = cell->vertex(l);
if (compute_jacobians)
- for (unsigned int point=0; point<n_points; ++point)
+ switch (dim)
{
- const double xi = unit_points[point](0);
- const double eta= unit_points[point](1);
+ case 1:
+ for (unsigned int point=0; point<n_points; ++point)
+ jacobians[point][0][0] = 1./(vertices[1](0)-vertices[0](0));
+ break;
+
+ case 2:
+ {
+ for (unsigned int point=0; point<n_points; ++point)
+ {
+ const double xi = unit_points[point](0);
+ const double eta= unit_points[point](1);
- const double t6 = vertices[0](0)*vertices[3](1);
- const double t8 = vertices[2](0)*xi;
- const double t10 = vertices[1](0)*eta;
- const double t12 = vertices[3](0)*vertices[1](1);
- const double t16 = vertices[3](0)*xi;
- const double t20 = vertices[0](0)*vertices[1](1);
- const double t22 = vertices[0](0)*vertices[2](1);
- const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
- t12*eta-vertices[3](0)*vertices[2](1)*eta-
- t16*vertices[0](1)+t16*vertices[1](1)-t12+
- vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
- const double t28 = vertices[1](0)*vertices[3](1);
- const double t31 = vertices[2](0)*eta;
- const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
- t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
- t31*vertices[3](1)+t20-t6-vertices[1](0)*
- vertices[0](1)+t28-t22*xi;
- const double t38 = 1/(t24+t36);
-
- jacobians[point][0][0] = (-vertices[0](1)+vertices[0](1)*xi-
- vertices[1](1)*xi+vertices[2](1)*xi+
- vertices[3](1)-vertices[3](1)*xi)*t38;
- jacobians[point][0][1] = -(-vertices[0](0)+vertices[0](0)*xi-
- vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
- jacobians[point][1][0] = -(-vertices[0](1)+vertices[0](1)*eta+
- vertices[1](1)-vertices[1](1)*eta+
- vertices[2](1)*eta-vertices[3](1)*eta)*t38;
- jacobians[point][1][1] = (-vertices[0](0)+vertices[0](0)*eta+
- vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
+ const double t6 = vertices[0](0)*vertices[3](1);
+ const double t8 = vertices[2](0)*xi;
+ const double t10 = vertices[1](0)*eta;
+ const double t12 = vertices[3](0)*vertices[1](1);
+ const double t16 = vertices[3](0)*xi;
+ const double t20 = vertices[0](0)*vertices[1](1);
+ const double t22 = vertices[0](0)*vertices[2](1);
+ const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
+ t12*eta-vertices[3](0)*vertices[2](1)*eta-
+ t16*vertices[0](1)+t16*vertices[1](1)-t12+
+ vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
+ const double t28 = vertices[1](0)*vertices[3](1);
+ const double t31 = vertices[2](0)*eta;
+ const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
+ t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
+ t31*vertices[3](1)+t20-t6-vertices[1](0)*
+ vertices[0](1)+t28-t22*xi;
+ const double t38 = 1/(t24+t36);
+
+ jacobians[point][0][0] = (-vertices[0](1)+vertices[0](1)*xi-
+ vertices[1](1)*xi+vertices[2](1)*xi+
+ vertices[3](1)-vertices[3](1)*xi)*t38;
+ jacobians[point][0][1] = -(-vertices[0](0)+vertices[0](0)*xi-
+ vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
+ jacobians[point][1][0] = -(-vertices[0](1)+vertices[0](1)*eta+
+ vertices[1](1)-vertices[1](1)*eta+
+ vertices[2](1)*eta-vertices[3](1)*eta)*t38;
+ jacobians[point][1][1] = (-vertices[0](0)+vertices[0](0)*eta+
+ vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
+ };
+
+ break;
+ };
+
+ default:
+ // not implemented at present
+ Assert (false, ExcNotImplemented());
};
};
template <int dim>
-void FELinear<dim>::get_support_points (const typename DoFHandler<dim>::cell_iterator &cell,
- const Boundary<dim> &,
- vector<Point<dim> > &support_points) const {
+void
+FELinear<dim>::get_support_points (const typename DoFHandler<dim>::cell_iterator &cell,
+ const Boundary<dim> &,
+ vector<Point<dim> > &support_points) const {
Assert (support_points.size() == total_dofs,
ExcWrongFieldDimension (support_points.size(), total_dofs));
template <int dim>
-void FELinear<dim>::get_face_support_points (const typename DoFHandler<dim>::face_iterator &face,
- const Boundary<dim> &,
- vector<Point<dim> > &support_points) const {
+void
+FELinear<dim>::get_face_support_points (const typename DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &,
+ vector<Point<dim> > &support_points) const {
Assert ((support_points.size() == dofs_per_face) &&
(support_points.size() == GeometryInfo<dim>::vertices_per_face),
ExcWrongFieldDimension (support_points.size(),
template <>
inline
-Point<1>
+Tensor<1,1>
FELinearMapping<1>::shape_grad_transform(const unsigned int i,
const Point<1>&) const
{
vector<Point<1> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<1> > > &shape_gradients_transform,
+ const vector<vector<Tensor<1,1> > > &shape_gradients_transform,
const Boundary<1> &boundary) const {
// simply pass down
FiniteElement<1>::fill_fe_values (cell, unit_points,
template <>
inline
-Point<2>
+Tensor<1,2>
FELinearMapping<2>::shape_grad_transform (const unsigned int i,
const Point<2>& p) const
{
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
- const vector<vector<Point<dim> > > &/*shape_grad_transform*/,
+ const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
is multiplied to the unit cell gradient *from the right*! be very careful
with these things.
- The following little program tests the correct behaviour (you have to find
- out the right include files, I tested it within a whole project with far
- more include files than necessary):
+ The following little program tests the correct behaviour:
-------------------------------------------
+ #include <grid/tria.h>
+ #include <grid/tria_boundary.h>
+ #include <grid/dof.h>
+ #include <fe/fe_values.h>
+ #include <fe/fe_lib.lagrange.h>
+ #include <base/quadrature_lib.h>
+ #include <grid/tria_iterator.h>
+ #include <grid/dof_accessor.h>
+ #include <lac/dvector.h>
+
int main () {
Triangulation<2> tria;
tria.create_hypercube (0,1);
*/
if (compute_jacobians)
- for (unsigned int point=0; point<n_points; ++point)
- switch (dim)
+ switch (dim)
+ {
+ case 1:
+ for (unsigned int point=0; point<n_points; ++point)
+ jacobians[point][0][0] = 1./(vertices[1](0)-vertices[0](0));
+ break;
+
+ case 2:
{
- case 2:
- {
-
- const double xi = unit_points[point](0);
- const double eta= unit_points[point](1);
+ for (unsigned int point=0; point<n_points; ++point)
+ {
+ const double xi = unit_points[point](0);
+ const double eta= unit_points[point](1);
- const double t6 = vertices[0](0)*vertices[3](1);
- const double t8 = vertices[2](0)*xi;
- const double t10 = vertices[1](0)*eta;
- const double t12 = vertices[3](0)*vertices[1](1);
- const double t16 = vertices[3](0)*xi;
- const double t20 = vertices[0](0)*vertices[1](1);
- const double t22 = vertices[0](0)*vertices[2](1);
- const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
- t12*eta-vertices[3](0)*vertices[2](1)*eta-
- t16*vertices[0](1)+t16*vertices[1](1)-t12+
- vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
- const double t28 = vertices[1](0)*vertices[3](1);
- const double t31 = vertices[2](0)*eta;
- const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
- t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
- t31*vertices[3](1)+t20-t6-vertices[1](0)*
- vertices[0](1)+t28-t22*xi;
- const double t38 = 1/(t24+t36);
-
- jacobians[point][0][0] = (-vertices[0](1)+vertices[0](1)*xi-
- vertices[1](1)*xi+vertices[2](1)*xi+
- vertices[3](1)-vertices[3](1)*xi)*t38;
- jacobians[point][0][1] = -(-vertices[0](0)+vertices[0](0)*xi-
- vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
- jacobians[point][1][0] = -(-vertices[0](1)+vertices[0](1)*eta+
- vertices[1](1)-vertices[1](1)*eta+
- vertices[2](1)*eta-vertices[3](1)*eta)*t38;
- jacobians[point][1][1] = (-vertices[0](0)+vertices[0](0)*eta+
- vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
-
- break;
- };
-
- default:
- // not implemented at present,
- // because of the changes above
- Assert (false, ExcNotImplemented());
+ const double t6 = vertices[0](0)*vertices[3](1);
+ const double t8 = vertices[2](0)*xi;
+ const double t10 = vertices[1](0)*eta;
+ const double t12 = vertices[3](0)*vertices[1](1);
+ const double t16 = vertices[3](0)*xi;
+ const double t20 = vertices[0](0)*vertices[1](1);
+ const double t22 = vertices[0](0)*vertices[2](1);
+ const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
+ t12*eta-vertices[3](0)*vertices[2](1)*eta-
+ t16*vertices[0](1)+t16*vertices[1](1)-t12+
+ vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
+ const double t28 = vertices[1](0)*vertices[3](1);
+ const double t31 = vertices[2](0)*eta;
+ const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
+ t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
+ t31*vertices[3](1)+t20-t6-vertices[1](0)*
+ vertices[0](1)+t28-t22*xi;
+ const double t38 = 1/(t24+t36);
+
+ jacobians[point][0][0] = (-vertices[0](1)+vertices[0](1)*xi-
+ vertices[1](1)*xi+vertices[2](1)*xi+
+ vertices[3](1)-vertices[3](1)*xi)*t38;
+ jacobians[point][0][1] = -(-vertices[0](0)+vertices[0](0)*xi-
+ vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
+ jacobians[point][1][0] = -(-vertices[0](1)+vertices[0](1)*eta+
+ vertices[1](1)-vertices[1](1)*eta+
+ vertices[2](1)*eta-vertices[3](1)*eta)*t38;
+ jacobians[point][1][1] = (-vertices[0](0)+vertices[0](0)*eta+
+ vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
+ };
+
+ break;
};
+
+ default:
+ // not implemented at present
+ Assert (false, ExcNotImplemented());
+ };
template <int dim>
-void FEValuesBase<dim>::get_function_grads (const dVector &fe_function,
- vector<Point<dim> > &gradients) const {
+void FEValuesBase<dim>::get_function_grads (const dVector &fe_function,
+ vector<Tensor<1,dim> > &gradients) const {
Assert (gradients.size() == n_quadrature_points,
ExcWrongVectorSize(gradients.size(), n_quadrature_points));
present_cell->get_dof_values (fe_function, dof_values);
// initialize with zero
- fill_n (gradients.begin(), n_quadrature_points, Point<dim>());
+ fill_n (gradients.begin(), n_quadrature_points, Tensor<1,dim>());
// add up contributions of trial
// functions
template <int dim>
-const Point<dim> & FEValuesBase<dim>::quadrature_point (const unsigned int i) const {
+const Point<dim> &
+FEValuesBase<dim>::quadrature_point (const unsigned int i) const {
Assert (i<quadrature_points.size(), ExcInvalidIndex(i,quadrature_points.size()));
Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
template <int dim>
-const Point<dim> & FEValuesBase<dim>::support_point (const unsigned int i) const {
+const Point<dim> &
+FEValuesBase<dim>::support_point (const unsigned int i) const {
Assert (i<support_points.size(), ExcInvalidIndex(i, support_points.size()));
Assert (update_flags & update_support_points, ExcAccessToUninitializedField());
unit_shape_2nd_derivatives(fe.total_dofs,
vector<Tensor<2,dim> >(quadrature.n_quadrature_points)),
unit_shape_gradients_transform(fe.n_transform_functions,
- vector<Point<dim> >(quadrature.n_quadrature_points)),
+ vector<Tensor<1,dim> >(quadrature.n_quadrature_points)),
unit_quadrature_points(quadrature.get_quad_points())
{
Assert ((update_flags & update_normal_vectors) == false,
vector<vector<Tensor<2,dim> > >(n_dofs,
vector<Tensor<2,dim> >(n_q_points))),
unit_shape_gradients_transform (n_faces_or_subfaces,
- vector<vector<Point<dim> > >(n_transform_functions,
- vector<Point<dim> >(n_q_points))),
+ vector<vector<Tensor<1,dim> > >(n_transform_functions,
+ vector<Tensor<1,dim> >(n_q_points))),
unit_face_quadrature_points (n_q_points, Point<dim-1>()),
unit_quadrature_points (n_faces_or_subfaces,
vector<Point<dim> >(n_q_points, Point<dim>())),
template <int dim>
-const Point<dim> & FEFaceValuesBase<dim>::normal_vector (const unsigned int i) const {
+const Point<dim> &
+FEFaceValuesBase<dim>::normal_vector (const unsigned int i) const {
Assert (i<normal_vectors.size(), ExcInvalidIndex(i, normal_vectors.size()));
Assert (update_flags & update_normal_vectors,
ExcAccessToUninitializedField());
{
fill_n (shape_gradients[i].begin(),
n_quadrature_points,
- Point<dim>());
+ Tensor<1,dim>());
for (unsigned int j=0; j<n_quadrature_points; ++j)
for (unsigned int s=0; s<dim; ++s)
// (grad psi)_s =
{
fill_n (shape_gradients[i].begin(),
n_quadrature_points,
- Point<dim>());
+ Tensor<1,dim>());
for (unsigned int j=0; j<n_quadrature_points; ++j)
for (unsigned int s=0; s<dim; ++s)
// (grad psi)_s =
//
// let psi be a short name for
// [a grad u_h]
- vector<Point<dim> > psi(n_q_points);
+ vector<Tensor<1,dim> > psi(n_q_points);
fe_face_values_cell.get_function_grads (solution, psi);
// get a list of the gradients of
// the finite element solution
// restricted to the neighbor cell
- vector<Point<dim> > neighbor_psi (n_q_points);
+ vector<Tensor<1,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> >());
+ minus<Tensor<1,dim> >());
};
//
// let psi be a short name for
// [a grad u_h]
- vector<Point<dim> > psi(n_q_points);
+ vector<Tensor<1,dim> > psi(n_q_points);
// store which number #cell# has in the
// list of neighbors of #neighbor#
// restrict the finite element on the
// neighbor cell to the common #subface#.
// store the gradient in #neighbor_psi#
- vector<Point<dim> > neighbor_psi (n_q_points);
+ vector<Tensor<1,dim> > neighbor_psi (n_q_points);
fe_face_values.reinit (neighbor_child, neighbor_neighbor, boundary);
fe_face_values.get_function_grads (solution, neighbor_psi);
transform (psi.begin(), psi.end(),
neighbor_psi.begin(),
psi.begin(),
- minus<Point<dim> >());
+ minus<Tensor<1,dim> >());
// next we have to multiply this with
// the normal vector. Since we have
template <int dim>
-inline double sqr_point (const Point<dim> &p) {
- return p.square();
+inline double sqr_point (const Tensor<1,dim> &p) {
+ return p * p;
};
// same procedure as above, but now
// psi is a vector of gradients
const unsigned int n_q_points = q.n_quadrature_points;
- vector<Point<dim> > psi (n_q_points);
+ vector<Tensor<1,dim> > psi (n_q_points);
// in praxi: first compute
// exact fe_function vector
// fe_function
if (true)
{
- vector<Point<dim> > function_grads (n_q_points, Point<dim>());
+ vector<Tensor<1,dim> > function_grads (n_q_points, Tensor<1,dim>());
fe_values.get_function_grads (fe_function, function_grads);
transform (psi.begin(), psi.end(),
function_grads.begin(),
psi.begin(),
- minus<Point<dim> >());
+ minus<Tensor<1,dim> >());
};
// take square of integrand
vector<double> psi_square (psi.size(), 0.0);