]> https://gitweb.dealii.org/ - dealii.git/commitdiff
More cell transform stuff.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 6 Mar 1998 13:54:50 +0000 (13:54 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 6 Mar 1998 13:54:50 +0000 (13:54 +0000)
git-svn-id: https://svn.dealii.org/trunk@33 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_lib.linear.cc

index 0c392f06084245ebe8070635f8ac93e8c5e049af..f1b7be3b181eeb1a2df4243a049c2be44027afd8 100644 (file)
@@ -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<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> ();
@@ -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; 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);
+    };
 };
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.