]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add support for the exact evaluation of mass matrices to the finite element classes.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 25 May 1998 14:39:32 +0000 (14:39 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 25 May 1998 14:39:32 +0000 (14:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@352 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.lagrange.h
deal.II/deal.II/source/fe/fe_lib.linear.cc

index 6e3fad1ca8edea5fb6df2791c163b1d09b0dda8a..523e5fc65f4a5c3f22154b393d816aa447c9f46b 100644 (file)
@@ -342,7 +342,7 @@ struct FiniteElementBase : public FiniteElementData<dim> {
  * implementation.
  *
  *
- * \subsection{Finite Elements in one dimension}
+ * \subsection{Finite elements in one dimension}
  *
  * Finite elements in one dimension need only set the #restriction# and
  * #prolongation# matrices in #FiniteElementBase<1>#. The constructor of
@@ -692,7 +692,7 @@ class FiniteElement : public FiniteElementBase<dim> {
                                     /**
                                      * Compute the off-points of the finite
                                      * element basis functions on the given
-                                     * cell.
+                                     * cell in real space.
                                      *
                                      * This function implements a subset of
                                      * the information delivered by the
@@ -835,6 +835,85 @@ class FiniteElement : public FiniteElementBase<dim> {
                                     const vector<Point<dim-1> > &unit_points,
                                     vector<Point<dim> >         &normal_vectors) const =0;
 
+                                    /**
+                                     * Fill in the given matrix with the local
+                                     * mass matrix. The mass matrix must be
+                                     * exactly computed, not using a
+                                     * quadrature, which may be done using
+                                     * an equation object and an assembler,
+                                     * as is done for the Laplace matrix
+                                     * in the #MatrixTools# class for example.
+                                     *
+                                     * The exact integration is possible since
+                                     * an exact representation for the Jacobi
+                                     * determinant exists in all known cases of
+                                     * iso- or subparametric mappings. For
+                                     * example, usually the point in real
+                                     * space $\vec x$ referring to the point
+                                     * $\vec \xi$ on the unit cell is given
+                                     * by $\vec x = \sum_i \vec p_i \phi_i(\vec \xi)$,
+                                     * where the sum is over all basis functions
+                                     * $\phi_i$ and $\vec p_i$ are the points
+                                     * in real space where the basis function
+                                     * $\phi_i$ is located. The Jacobi
+                                     * determinant is the given by
+                                     * $|det J| = |\frac{\partial\vec x}{\partial\vec\xi}$,
+                                     * which can be evaluated in closed form.
+                                     * The mass matrix then is given by
+                                     * $m_{ij} = \int_{\hat K} \phi_i(\vec\xi)
+                                     * \phi_j(\vec\xi) |det J| d\xi$, where
+                                     * $\hat K$ is the unit cell. The integrand
+                                     * obviously is a polynom and can thus
+                                     * easily be integrated analytically, so
+                                     * the computation of the local mass matrix
+                                     * is reduced to the computation of a
+                                     * weighted evaluation of a polynom in
+                                     * the coordinates of the ansatz points
+                                     * in real space (for linear mappings,
+                                     * these are the corner points, for
+                                     * quadratic mappings also the center of
+                                     * mass and the edge and face centers).
+                                     * For example, in one space dimension,
+                                     * the Jacobi determinant simply is $h$,
+                                     * the size of the cell, and the integral
+                                     * over the two basis functions can easily
+                                     * be calculated with a pen and a sheet of
+                                     * paper. The actual computation on this
+                                     * matrix then is simply a scaling of a
+                                     * known and constant matrix by $h$.
+                                     *
+                                     * The functions which override this one
+                                     * may make assumptions on the sign of
+                                     * the determinant if stated in the
+                                     * documentation, but should check for
+                                     * them in debug mode. For that purpose,
+                                     * an exception with the longish name
+                                     * #ExcJacobiDeterminantHasWrongSign#
+                                     * is declared.
+                                     *
+                                     * The function takes a #DoFHandler#
+                                     * iterator, which provides a superset
+                                     * of information to the geometrical
+                                     * information needed for the computations.
+                                     * The additional data should not be
+                                     * used, however a #DoFHandler# iterator
+                                     * was preferred over a #Triangulation#
+                                     * iterator since this is what usually
+                                     * is available in places where this
+                                     * function is called.
+                                     *
+                                     * The cell matrix is assumed to be of
+                                     * the right size already. Functions
+                                     * of derived classes shall be implemented
+                                     * in a way as to overwrite the previous
+                                     * contents of the matrix, so it need not
+                                     * be necessary to clear the matrix before
+                                     * use with this function.
+                                     */
+    virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+                                       const Boundary<dim> &boundary, 
+                                       dFMatrix            &local_mass_matrix) const =0;
+    
                                     /**
                                      * Exception
                                      */
@@ -847,6 +926,10 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      * Exception
                                      */
     DeclException0 (ExcBoundaryFaceUsed);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcJacobiDeterminantHasWrongSign);
 };
 
 
index 93dc2a28916d04cd01cd2badc93d895d9fe30cd9..acb5c59cbd24478cd6d1bc4f770c75a48c81a607 100644 (file)
@@ -148,6 +148,14 @@ class FELinear : public FiniteElement<dim> {
                                     const unsigned int           subface_no,
                                     const vector<Point<dim-1> > &unit_points,
                                     vector<Point<dim> >         &normal_vectors) const;    
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+                                       const Boundary<dim> &boundary,
+                                       dFMatrix &local_mass_matrix) const;
 };
 
 
@@ -251,6 +259,14 @@ class FEQuadratic : public FiniteElement<dim> {
                                     const unsigned int           face_no,
                                     const vector<Point<dim-1> > &unit_points,
                                     vector<Point<dim> >         &normal_vectors) const;    
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+                                       const Boundary<dim> &boundary,
+                                       dFMatrix &local_mass_matrix) const;
 };
 
 
@@ -354,6 +370,14 @@ class FECubic : public FiniteElement<dim> {
                                     const unsigned int           face_no,
                                     const vector<Point<dim-1> > &unit_points,
                                     vector<Point<dim> >         &normal_vectors) const;    
+
+                                    /**
+                                     * Refer to the base class for detailed
+                                     * information on this function.
+                                     */
+    virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
+                                       const Boundary<dim> &boundary,
+                                       dFMatrix &local_mass_matrix) const;
 };
 
 
index 22d559484e7621a3935ea14c9cef7bdbd4f99acb..aeb5c4707e0fdd2ffea12b32f35c84ec78add996 100644 (file)
@@ -154,6 +154,22 @@ void FELinear<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
 
 
 
+template <>
+void FELinear<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell,
+                                        const Boundary<1> &,
+                                        dFMatrix &local_mass_matrix) const {
+  Assert (local_mass_matrix.n() == total_dofs,
+         ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs));
+  Assert (local_mass_matrix.m() == total_dofs,
+         ExcWrongFieldDimension(local_mass_matrix.m(),total_dofs));
+
+  const double h = cell->vertex(1)(0) - cell->vertex(0)(0);
+  Assert (h>0, ExcJacobiDeterminantHasWrongSign());
+
+  local_mass_matrix(0,0) = local_mass_matrix(1,1) = 1./3.*h;
+  local_mass_matrix(0,1) = local_mass_matrix(1,0) = 1./6.*h;
+};
+
 
 
 template <>
@@ -260,7 +276,7 @@ FELinear<2>::shape_value (const unsigned int i,
     case 3: return (1.-p(0)) * p(1);
     }
   return 0.;
-}
+};
 
 
 
@@ -278,7 +294,111 @@ FELinear<2>::shape_grad (const unsigned int i,
     case 3: return Point<2> (-p(1), 1.-p(0));
     }
   return Point<2> ();
-}
+};
+
+
+
+template <>
+void FELinear<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell,
+                                        const Boundary<2> &,
+                                        dFMatrix &local_mass_matrix) const {
+  Assert (local_mass_matrix.n() == total_dofs,
+         ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs));
+  Assert (local_mass_matrix.m() == total_dofs,
+         ExcWrongFieldDimension(local_mass_matrix.m(),total_dofs));
+
+/* Get the computation of the local mass matrix by these lines in maple:
+
+   x_real := sum(x[i]*phi[i], i=0..3);
+   y_real := sum(y[i]*phi[i], i=0..3);
+   phi[0] := (1-xi)*(1-eta);
+   phi[1] := xi*(1-eta);
+   phi[2] := xi*eta;
+   phi[3] := (1-xi)*eta;
+   detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi);
+
+   m := proc (i,j)  int( int(phi[i]*phi[j]*detJ, xi=0..1), eta=0..1); end;
+
+   M := array(0..3,0..3);
+   for i from 0 to 3 do
+     for j from 0 to 3 do
+       M[i,j] := m(i,j);
+     od;
+   od;
+
+   readlib(C);
+   C(M, optimized);
+*/
+
+  const double x[4] = { cell->vertex(0)(0),
+                       cell->vertex(1)(0),
+                       cell->vertex(2)(0),
+                       cell->vertex(3)(0)  };
+  const double y[4] = { cell->vertex(0)(1),
+                       cell->vertex(1)(1),
+                       cell->vertex(2)(1),
+                       cell->vertex(3)(1)  };
+  
+/* check that the Jacobi determinant
+
+    t0 = (-x[0]*(1.0-eta)+x[1]*(1.0-eta)+x[2]*eta-x[3]*eta) *
+         (-y[0]*(1.0-xi)-y[1]*xi+y[2]*xi+y[3]*(1.0-xi))        -
+        (-x[0]*(1.0-xi)-x[1]*xi+x[2]*xi+x[3]*(1.0-xi)) *
+        (-y[0]*(1.0-eta)+y[1]*(1.0-eta)+y[2]*eta-y[3]*eta)
+
+   has the right sign.  
+        
+   We do not attempt to check its (hopefully) positive sign at all points
+   on the unit cell, but we check that it is positive in the four corners,
+   which is sufficient since $det J$ is a bilinear function.
+*/
+  Assert ((-x[0]+x[1])*(-y[0]+y[3])-(-x[0]+x[3])*(-y[0]+y[1]),  // xi=eta=0
+         ExcJacobiDeterminantHasWrongSign());
+  Assert ((x[2]-x[3])*(-y[0]+y[3])-(-x[0]+x[3])*(y[2]-y[3]),    // xi=0, eta=1
+         ExcJacobiDeterminantHasWrongSign());
+  Assert ((x[2]-x[3])*(-y[1]+y[2])-(-x[1]+x[2])*(y[2]-y[3]),    // xi=eta=1
+         ExcJacobiDeterminantHasWrongSign());
+  Assert ((-x[0]+x[1])*(-y[1]+y[2])-(-x[1]+x[2])*(-y[0]+y[1]),  // xi=1, eta=0
+         ExcJacobiDeterminantHasWrongSign());
+  
+  const double t1 = x[1]*y[3],
+              t2 = x[1]*y[2],
+              t3 = x[1]*y[0],
+              t4 = x[0]*y[3],
+              t5 = x[0]*y[1],
+              t6 = x[2]*y[3],
+              t7 = x[3]*y[0],
+              t8 = x[2]*y[1],
+              t9 = x[3]*y[2],
+             t10 = x[3]*y[1],
+             t12 = x[0]*y[2],
+             t13 = x[2]*y[0],
+             t14 = t1/72+t2/36-t3/24-t4/36-t12/72+t5/24+t6/72
+                   +t7/36-t8/36-t9/72-t10/72+t13/72,
+             t15 = t2/72-t3/72-t4/72+t5/72+t6/72+t7/72-t8/72-t9/72,
+             t16 = t1/72+t2/72-t3/36-t4/24+t12/72+t5/36+t6/36
+                   +t7/24-t8/72-t9/36-t10/72-t13/72,
+             t18 = -t1/72+t2/24-t3/36-t4/72-t12/72+t5/36+t6/36
+                   +t7/72-t8/24-t9/36+t10/72+t13/72,
+             t20 = -t1/72+t12/72+t2/36+t5/72-t3/72+t6/24
+                   -t9/24-t13/72+t10/72-t4/36+t7/36-t8/36;
+  local_mass_matrix(0,0) = t1/18+t2/36-t3/12-t4/12+t5/12+t6/36+t7/12-t8/36-t9/36-t10/18;
+  local_mass_matrix(0,1) = t14;
+  local_mass_matrix(0,2) = t15;
+  local_mass_matrix(0,3) = t16;
+  local_mass_matrix(1,0) = t14;
+  local_mass_matrix(1,1) = t2/12-t3/12-t4/36-t12/18+t5/12+t6/36+t7/36-t8/12-t9/36+t13/18;
+  local_mass_matrix(1,2) = t18;
+  local_mass_matrix(1,3) = t15;
+  local_mass_matrix(2,0) = t15;
+  local_mass_matrix(2,1) = t18;
+  local_mass_matrix(2,2) = -t1/18+t2/12+t5/36-t3/36+t6/12-t9/12+t10/18-t4/36+t7/36-t8/12;
+  local_mass_matrix(2,3) = t20;
+  local_mass_matrix(3,0) = t16;
+  local_mass_matrix(3,1) = t15;
+  local_mass_matrix(3,2) = t20;
+  local_mass_matrix(3,3) = t12/18+t2/36+t5/36-t3/36+t6/12-t9/12-t13/18-t4/12+t7/12-t8/36;
+};
 
 
 
@@ -599,6 +719,7 @@ void FEQuadratic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
 
 
 
+
 template <>
 FEQuadratic<2>::FEQuadratic () :
                FiniteElement<2> (1, 1, 1)
@@ -729,6 +850,17 @@ void FEQuadratic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator
 
 
 
+template <int dim>
+void FEQuadratic<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &,
+                                             const Boundary<dim> &,
+                                             dFMatrix &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
+
 
 
 
@@ -936,6 +1068,16 @@ void FECubic<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cel
 
 
 
+template <int dim>
+void FECubic<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &,
+                                         const Boundary<dim> &,
+                                         dFMatrix &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
 // explicit instantiations
 
 template class FELinear<1>;

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.