From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 6 Mar 1998 10:52:31 +0000 (+0000)
Subject: Cell transform stuff.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cd2debe824006b1cbb5ce031b38071d0b67aca37;p=dealii-svn.git

Cell transform stuff.


git-svn-id: https://svn.dealii.org/trunk@29 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc
index 501a17d0b4..79db61b61c 100644
--- a/deal.II/deal.II/source/fe/fe.cc
+++ b/deal.II/deal.II/source/fe/fe.cc
@@ -2,6 +2,10 @@
 
 #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
 
 
 
@@ -89,8 +93,11 @@ void FEValues<1>::reinit (const Triangulation<1>::cell_iterator &cell,
   
 				   // 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();
 };
 
 
@@ -120,10 +127,11 @@ void FEValues<2>::reinit (const Triangulation<2>::cell_iterator &cell,
 	      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();
 };
 
 
@@ -226,6 +234,23 @@ bool FiniteElement<1>::operator == (const FiniteElement<1> &f) const {
 
 
 
+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) &&
@@ -235,6 +260,15 @@ bool FiniteElement<2>::operator == (const FiniteElement<2> &f) const {
 
 
 
+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 -------------*/
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 a56521b238..0c392f0608 100644
--- a/deal.II/deal.II/source/fe/fe_lib.linear.cc
+++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc
@@ -1,7 +1,11 @@
 /* $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 () :
@@ -75,6 +79,18 @@ FELinear<1>::shape_grad(const unsigned int i,
 
 
 
+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)
 {
@@ -176,8 +192,8 @@ FELinear<2>::FELinear () :
 
 
 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)
@@ -193,8 +209,8 @@ FELinear<2>::shape_value(const unsigned int 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)
@@ -209,10 +225,56 @@ 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
+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)
 {
@@ -229,10 +291,23 @@ FEQuadratic<2>::FEQuadratic () :
 
 
 
+
 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) {};