#include <fe/fe.h>
#include <fe/quadrature.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+
+extern TriaIterator<1,CellAccessor<1> > __dummy2687; // for gcc2.8
// compute Jacobi determinants in
// quadrature points.
+ // refer to the general doc for
+ // why we take the inverse of the
+ // determinant
for (unsigned int i=0; i<quadrature_points.size(); ++i)
- JxW_values[i] = jacobi_matrices[i].determinant() * weights[i];
+ JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
};
unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
};
- // compute Jacobi determinants in
- // quadrature points.
+ // refer to the general doc for
+ // why we take the inverse of the
+ // determinant
for (unsigned int i=0; i<quadrature_points.size(); ++i)
- JxW_values[i] = jacobi_matrices[i].determinant() * weights[i];
+ JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
};
+void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
+ const vector<Point<1> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ vector<Point<1> > &points) const {
+ // local mesh width
+ double h=(cell->vertex(1)(0) - cell->vertex(0)(0));
+
+ unsigned int n_points = unit_points.size();
+ for (unsigned int i=0; i<n_points; ++i)
+ {
+ jacobians[i](0,0) = 1./h;
+ points[i] = cell->vertex(0) + h*unit_points[i];
+ };
+};
+
+
+
bool FiniteElement<2>::operator == (const FiniteElement<2> &f) const {
return ((dofs_per_vertex == f.dofs_per_vertex) &&
(dofs_per_line == f.dofs_per_line) &&
+void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
+ const vector<Point<2> > &,
+ vector<dFMatrix> &,
+ vector<Point<2> > &) const {
+ Assert (false, ExcPureFunctionCalled());
+};
+
+
+
/*------------------------------- Explicit Instantiations -------------*/
/* $Id$ */
#include <fe/fe_lib.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+extern TriaIterator<1,CellAccessor<1> > __dummy2687; // for gcc2.8
+extern TriaIterator<2,CellAccessor<2> > __dummy2688; // for gcc2.8
FELinear<1>::FELinear () :
+void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
+ const vector<Point<1> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ vector<Point<1> > &points) const {
+ // simply pass down
+ FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+};
+
+
+
+
+
FELinear<2>::FELinear () :
FiniteElement<2> (1, 0, 0)
{
double
-FELinear<2>::shape_value(const unsigned int i,
- const Point<2>& p) const
+FELinear<2>::shape_value (const unsigned int i,
+ const Point<2>& p) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
switch (i)
Point<2>
-FELinear<2>::shape_grad(const unsigned int i,
- const Point<2>& p) const
+FELinear<2>::shape_grad (const unsigned int i,
+ const Point<2>& p) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
switch (i)
+// this function may be generalised to three or more dimensions with gcc2.8
+// you will have to change th number of vertices
+void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
+ const vector<Point<2> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ vector<Point<2> > &points) const {
+ const unsigned int dim=2;
+ const unsigned int n_vertices=4;
+
+ unsigned int n_points=unit_points.size();
+
+ // initialize points to zero
+ for (unsigned int i=0; i<n_points; ++i)
+ points[i] = Point<dim> ();
+
+ // note: let x_l be the vector of the
+ // lth quadrature point in real space and
+ // xi_l that on the unit cell, let further
+ // p_j be the vector of the jth vertex
+ // of the cell in real space and
+ // N_j(xi_l) be the value of the associated
+ // basis function at xi_l, then
+ // x_l(xi_l) = sum_j p_j N_j(xi_l)
+ for (unsigned int j=0; j<n_vertices; ++j)
+ for (unsigned int l=0; l<n_points; ++l)
+ points[l] += cell->vertex(j) * shape_value(j, unit_points[l]);
+
+// computation of jacobian still missing
+};
+
+
+
+
+
+
FEQuadratic<1>::FEQuadratic () :
FiniteElement<1> (1, 1) {};
+
+void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
+ const vector<Point<1> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ vector<Point<1> > &points) const {
+ // simply pass down
+ FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+};
+
+
+
FEQuadratic<2>::FEQuadratic () :
FiniteElement<2> (1, 1, 1)
{
+
FECubic<1>::FECubic () :
FiniteElement<1> (1, 2) {};
+
+void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
+ const vector<Point<1> > &unit_points,
+ vector<dFMatrix> &jacobians,
+ vector<Point<1> > &points) const {
+ // simply pass down
+ FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+};
+
+
+
+
FECubic<2>::FECubic () :
FiniteElement<2> (1, 2, 4) {};