From 73aa7bb3cfeb90eb3d97fa5f6811efe45eac8fc0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 6 Mar 1998 13:54:50 +0000 Subject: [PATCH] More cell transform stuff. git-svn-id: https://svn.dealii.org/trunk@33 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_lib.linear.cc | 52 +++++++++++++++++++--- 1 file changed, 47 insertions(+), 5 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index 0c392f0608..f1b7be3b18 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -226,16 +226,22 @@ FELinear<2>::shape_grad (const unsigned int i, // 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 > &unit_points, vector &jacobians, vector > &points) const { const unsigned int dim=2; const unsigned int n_vertices=4; - unsigned int n_points=unit_points.size(); + // recomputation of values + Point vertices[n_vertices]; + for (unsigned int l=0; lvertex(l); + + // initialize points to zero for (unsigned int i=0; i (); @@ -250,9 +256,45 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, // x_l(xi_l) = sum_j p_j N_j(xi_l) for (unsigned int j=0; jvertex(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 gradient + = FELinear::shape_grad (s, unit_points[l]); + for (unsigned int i=0; i