From: Wolfgang Bangerth Date: Mon, 25 May 1998 14:39:32 +0000 (+0000) Subject: Add support for the exact evaluation of mass matrices to the finite element classes. X-Git-Tag: v8.0.0~22902 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=29eb5e27822d3a5382dba0bbf2ebc29e2f6655a6;p=dealii.git Add support for the exact evaluation of mass matrices to the finite element classes. git-svn-id: https://svn.dealii.org/trunk@352 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 6e3fad1ca8..523e5fc65f 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -342,7 +342,7 @@ struct FiniteElementBase : public FiniteElementData { * 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 { /** * 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 { const vector > &unit_points, vector > &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::cell_iterator &cell, + const Boundary &boundary, + dFMatrix &local_mass_matrix) const =0; + /** * Exception */ @@ -847,6 +926,10 @@ class FiniteElement : public FiniteElementBase { * Exception */ DeclException0 (ExcBoundaryFaceUsed); + /** + * Exception + */ + DeclException0 (ExcJacobiDeterminantHasWrongSign); }; diff --git a/deal.II/deal.II/include/fe/fe_lib.lagrange.h b/deal.II/deal.II/include/fe/fe_lib.lagrange.h index 93dc2a2891..acb5c59cbd 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -148,6 +148,14 @@ class FELinear : public FiniteElement { const unsigned int subface_no, const vector > &unit_points, vector > &normal_vectors) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, + const Boundary &boundary, + dFMatrix &local_mass_matrix) const; }; @@ -251,6 +259,14 @@ class FEQuadratic : public FiniteElement { const unsigned int face_no, const vector > &unit_points, vector > &normal_vectors) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, + const Boundary &boundary, + dFMatrix &local_mass_matrix) const; }; @@ -354,6 +370,14 @@ class FECubic : public FiniteElement { const unsigned int face_no, const vector > &unit_points, vector > &normal_vectors) const; + + /** + * Refer to the base class for detailed + * information on this function. + */ + virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, + const Boundary &boundary, + dFMatrix &local_mass_matrix) const; }; 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 22d559484e..aeb5c4707e 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -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::get_normal_vectors (const DoFHandler::cell_iterator +template +void FEQuadratic::get_local_mass_matrix (const DoFHandler::cell_iterator &, + const Boundary &, + dFMatrix &) const { + Assert (false, ExcNotImplemented()); +}; + + + + + @@ -936,6 +1068,16 @@ void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &cel +template +void FECubic::get_local_mass_matrix (const DoFHandler::cell_iterator &, + const Boundary &, + dFMatrix &) const { + Assert (false, ExcNotImplemented()); +}; + + + + // explicit instantiations template class FELinear<1>;