* 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
/**
* 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
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
*/
* Exception
*/
DeclException0 (ExcBoundaryFaceUsed);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcJacobiDeterminantHasWrongSign);
};
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;
};
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;
};
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;
};
+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 <>
case 3: return (1.-p(0)) * p(1);
}
return 0.;
-}
+};
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;
+};
+
template <>
FEQuadratic<2>::FEQuadratic () :
FiniteElement<2> (1, 1, 1)
+template <int dim>
+void FEQuadratic<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &,
+ const Boundary<dim> &,
+ dFMatrix &) const {
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+
+
+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>;