// this function may be generalised to three or more dimensions with gcc2.8
-// you will have to change th number of vertices
+// you will have to change the number of vertices and the setting up of the
+// jacobi matrices
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();
+ // recomputation of values
+ Point<dim> vertices[n_vertices];
+ for (unsigned int l=0; l<n_vertices; ++l)
+ vertices[l] = cell->vertex(l);
+
+
// initialize points to zero
for (unsigned int i=0; i<n_points; ++i)
points[i] = Point<dim> ();
// 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
+ points[l] += vertices[j] * shape_value(j, unit_points[l]);
+
+/* jacobi matrices: compute d(x)/d(xi) and invert this
+ Let M(l) be the inverse of J at the quadrature point l, then
+ M_{ij}(l) = sum_s p_i(s) d(N_s(l))/d(xi_j)
+ where p_i(s) is the i-th coordinate of the s-th vertex vector,
+ N_s(l) is the value of the s-th vertex shape function at the
+ quadrature point l.
+
+ We could therefore write:
+ l=0..n_points-1
+ i=0..dim-1
+ j=0..dim-1
+ M_{ij}(l) = 0
+ s=0..n_vertices
+ M_{ij}(l) += p_i(s) d(N_s(l))/d(xi_j)
+
+ However, we rewrite the loops to only compute the gradient once for
+ each integration point and basis function.
+*/
+ dFMatrix M(dim,dim);
+ for (unsigned int l=0; l<n_points; ++l)
+ {
+ M.clear ();
+ for (unsigned int s=0; s<n_vertices; ++s)
+ {
+ // we want a linear transform and
+ // if we prepend the class name in
+ // front of the #shape_grad#, we
+ // need not use virtual function
+ // calls.
+ const Point<dim> gradient
+ = FELinear<dim>::shape_grad (s, unit_points[l]);
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ M(i,j) += vertices[s](i) * gradient(j);
+ };
+ jacobians[l].invert(M);
+ };
};