/**
* Store a copy of the finite element given
- * latest for the distribution of dofs.
+ * latest for the distribution of dofs. In
+ * fact, since the FE class itself has not
+ * much functionality, this object only
+ * stores numbers such as the number of
+ * dofs per line etc. Calling any of the
+ * more specific functions will result in
+ * an error about calling a pure virtual
+ * function.
*/
FiniteElement<dim> *selected_fe;
the real cell at the quadrature points and so on.
The Jacobian matrix is defined to be
- $$ J_{ij} = {d\xi_i \over d\x_j} $$
+ $$ J_{ij} = {d\xi_i \over dx_j} $$
which is the form needed to compute the gradient on the real cell from
the gradient on the unit cell. If we want to transform the area element
$dx dy$ from the real to the unit cell, we have to take the determinant of
* should really be pure, but then we could
* not make copies of a finite element
* object even if we did not intend to use
- * this function. Therefore, we ommit the
+ * this function. Therefore, we omit the
* #=0# signature and implement this function
* by throwing an exception.
* #p# is a point on the reference element.
* should really be pure, but then we could
* not make copies of a finite element
* object even if we did not intend to use
- * this function. Therefore, we ommit the
+ * this function. Therefore, we omit the
* #=0# signature and implement this function
* by throwing an exception.
* #p# is a point on the reference element,
* higher dimensions, it depends on the
* present fe and needs reimplementation
* by the user.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
*/
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
const bool compute_q_points) const;
+
+ /**
+ * Return the ansatz points this FE has
+ * on a face if a cell would have the
+ * given face as a side. This function is
+ * needed for higher order elements, if
+ * we want to use curved boundary
+ * approximations.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
+ */
+ virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
+ vector<Point<dim> > &ansatz_points) const;
/**
* Comparison operator. We also check for
* Exception
*/
DeclException0 (ExcNotImplemented);
-
+ /**
+ * Exception
+ */
+ DeclException2 (ExcWrongFieldDimension,
+ int, int,
+ << "The field has not the assumed dimension " << arg2
+ << ", but has " << arg1 << " elements.");
protected:
/**
* Have #N=2^dim# matrices keeping the
elements which are no more valid.
Consequence: make sure the copy constructor is correct.
+
+ This class should really be a pure one, with functions having the #=0#
+ signature. However, some instances of this class need to be floating around
+ anyhow (e.g. the #DoFHandler# class keeps a copy, only to have the values
+ of #dof_per*# available), so we do not make it pure but rather implement
+ those functions which should in fact be pure to throw an error.
*/
class FiniteElement<1> : public FiniteElementBase<1> {
public:
* distance to the origin. The standard
* implementation distributes the dofs on
* the line equidistantly.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
*/
virtual void fill_fe_values (const Triangulation<1>::cell_iterator &cell,
const vector<Point<1> > &unit_points,
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
const bool compute_q_points) const;
+
+ /**
+ * Return the ansatz points this FE has
+ * on a face if a cell would have the
+ * given face as a side. This function is
+ * needed for higher order elements, if
+ * we want to use curved boundary
+ * approximations.
+ *
+ * Question: is this function useful in 1D?
+ * At present it is not implemented.
+ */
+ virtual void face_ansatz_points (const Triangulation<1>::face_iterator &face,
+ vector<Point<1> > &ansatz_points) const;
};
If you want to extend this class (not by derivation, but by adding new
elements), see \Ref{FiniteElement<1>}
+
+ This class should really be a pure one, with functions having the #=0#
+ signature. However, some instances of this class need to be floating around
+ anyhow (e.g. the #DoFHandler# class keeps a copy, only to have the values
+ of #dof_per*# available), so we do not make it pure but rather implement
+ those functions which should in fact be pure to throw an error.
*/
class FiniteElement<2> : public FiniteElementBase<2> {
public:
* function is therefore not implemented
* by the FE<2> base class, but is made
* pure virtual.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
*/
virtual void fill_fe_values (const Triangulation<2>::cell_iterator &cell,
const vector<Point<2> > &unit_points,
const bool compute_ansatz_points,
vector<Point<2> > &q_points,
const bool compute_q_points) const;
+
+ /**
+ * Return the ansatz points this FE has
+ * on a face if a cell would have the
+ * given face as a side. This function is
+ * needed for higher order elements, if
+ * we want to use curved boundary
+ * approximations.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
+ */
+ virtual void face_ansatz_points (const Triangulation<2>::face_iterator &face,
+ vector<Point<2> > &ansatz_points) const;
};
* function is therefore not implemented
* by the FE<2> base class, but is made
* pure virtual.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
*/
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
const bool compute_q_points) const;
+
+ /**
+ * Return the ansatz points this FE has
+ * on a face if a cell would have the
+ * given face as a side. Since we have no
+ * degrees of freedom on the faces for
+ * the linear ansatz, the ansatz points are
+ * simply the vertices of the face.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
+ */
+ virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
+ vector<Point<dim> > &ansatz_points) const;
};
* function is therefore not implemented
* by the FE<2> base class, but is made
* pure virtual.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
*/
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
* function is therefore not implemented
* by the FE<2> base class, but is made
* pure virtual.
+ *
+ * The function assumes that the fields
+ * already have the right number of
+ * elements.
*/
virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
const vector<Point<dim> > &unit_points,
* See the general doc for more
* information.
*/
- virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc,
- map<int,double> &boundary_values) const;
+ virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc,
+ const FiniteElement<dim> &fe,
+ map<int,double> &boundary_values) const;
/**
* Exception
void apply_dirichlet_bc (dSMatrix &matrix,
dVector &solution,
dVector &right_hand_side,
- const DirichletBC &dirichlet_bc);
+ const DirichletBC &dirichlet_bc,
+ const FiniteElement<dim> &fe);
friend class Assembler<dim>;
};
JxW_values(quadrature.n_quadrature_points, 0),
quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
unit_quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
+ ansatz_points (fe.total_dofs, Point<dim>()),
jacobi_matrices (quadrature.n_quadrature_points,
dFMatrix(dim,dim)),
update_flags (update_flags)
+template <int dim>
+void FiniteElementBase<dim>::face_ansatz_points (const typename Triangulation<dim>::face_iterator &,
+ vector<Point<dim> > &) const {
+ Assert (false, ExcPureFunctionCalled());
+};
+
+
+
/*------------------------------- FiniteElement ----------------------*/
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
const bool compute_q_points) const {
+ Assert (jacobians.size() == unit_points.size(),
+ ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+ Assert (q_points.size() == unit_points.size(),
+ ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+ Assert (ansatz_points.size() == total_dofs,
+ ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
+
// local mesh width
const double h=(cell->vertex(1)(0) - cell->vertex(0)(0));
+void FiniteElement<1>::face_ansatz_points (const typename Triangulation<1>::face_iterator &,
+ vector<Point<1> > &) const {
+ // is this function useful in 1D?
+ Assert (false, ExcPureFunctionCalled());
+};
+
+
bool FiniteElement<2>::operator == (const FiniteElement<2> &f) const {
return ((dofs_per_vertex == f.dofs_per_vertex) &&
(dofs_per_line == f.dofs_per_line) &&
+void FiniteElement<2>::face_ansatz_points (const typename Triangulation<2>::face_iterator &,
+ vector<Point<2> > &) const {
+ // is this function useful in 1D?
+ Assert (false, ExcPureFunctionCalled());
+};
+
+
+
/*------------------------------- Explicit Instantiations -------------*/
#include <grid/tria_accessor.h>
+
+
FELinear<1>::FELinear () :
FiniteElement<1> (1, 0)
{
+void FELinear<1>::face_ansatz_points (const Triangulation<1>::face_iterator &,
+ vector<Point<1> > &) const {
+ Assert (false, ExcPureFunctionCalled());
+};
+
+
FELinear<2>::FELinear () :
const bool compute_ansatz_points,
vector<Point<2> > &q_points,
const bool compute_q_points) const {
+ Assert (jacobians.size() == unit_points.size(),
+ ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+ Assert (q_points.size() == unit_points.size(),
+ ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+ Assert (ansatz_points.size() == total_dofs,
+ ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
const unsigned int dim=2;
- const unsigned int n_vertices=4;
+ const unsigned int n_vertices=(1<<dim);
unsigned int n_points=unit_points.size();
Point<dim> vertices[n_vertices];
+template <int dim>
+void FELinear<dim>::face_ansatz_points (const typename Triangulation<dim>::face_iterator &face,
+ vector<Point<dim> > &ansatz_points) const {
+ Assert (ansatz_points.size() == (1<<(dim-1)),
+ typename FiniteElementBase<dim>::ExcWrongFieldDimension (ansatz_points.size(),
+ 1<<(dim-1)));
+
+ for (unsigned int vertex=0; vertex<(1<<(dim-1)); ++vertex)
+ ansatz_points[vertex] = face->vertex(vertex);
+};
+
+
+
+
void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
- const vector<Point<2> > &,
- vector<dFMatrix> &,
+ const vector<Point<2> > &unit_points,
+ vector<dFMatrix> &jacobians,
const bool,
- vector<Point<2> > &,
+ vector<Point<2> > &ansatz_points,
const bool,
- vector<Point<2> > &,
+ vector<Point<2> > &q_points,
const bool) const {
+ Assert (jacobians.size() == unit_points.size(),
+ ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+ Assert (q_points.size() == unit_points.size(),
+ ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+ Assert (ansatz_points.size() == total_dofs,
+ ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
- const vector<Point<2> > &,
- vector<dFMatrix> &,
+ const vector<Point<2> > &unit_points,
+ vector<dFMatrix> &jacobians,
const bool,
- vector<Point<2> > &,
+ vector<Point<2> > &ansatz_points,
const bool,
- vector<Point<2> > &,
+ vector<Point<2> > &q_points,
const bool) const {
+ Assert (jacobians.size() == unit_points.size(),
+ ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
+ Assert (q_points.size() == unit_points.size(),
+ ExcWrongFieldDimension(q_points.size(), unit_points.size()));
+ Assert (ansatz_points.size() == total_dofs,
+ ExcWrongFieldDimension(ansatz_points.size(), total_dofs));
+
Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
// apply Dirichlet bc as described
// in the docs
apply_dirichlet_bc (system_matrix, solution,
- right_hand_side, dirichlet_bc);
+ right_hand_side,
+ dirichlet_bc, fe);
// condense system matrix in-place
constraints.condense (system_matrix);
void ProblemBase<dim>::apply_dirichlet_bc (dSMatrix &matrix,
dVector &solution,
dVector &right_hand_side,
- const DirichletBC &dirichlet_bc) {
+ const DirichletBC &dirichlet_bc,
+ const FiniteElement<dim> &fe) {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
ExcInvalidBoundaryIndicator());
-
+
// first make up a list of dofs subject
// to any boundary condition and which
// value they take; if a node occurs
// the lines in 2D being subject to
// different bc's), the last value is taken
map<int,double> boundary_values;
- make_boundary_value_list (dirichlet_bc, boundary_values);
+ make_boundary_value_list (dirichlet_bc, fe, boundary_values);
map<int,double>::const_iterator dof, endd;
const unsigned int n_dofs = (unsigned int)matrix.m();
// zero, take the first main
// diagonal entry we can find, or
// one if no nonzero main diagonal
- // element exists.
+ // element exists. Normally, however,
+ // the main diagonal entry should
+ // not be zero.
+ //
+ // store the new rhs entry to make
+ // the gauss step more efficient
+ double new_rhs;
if (matrix.diag_element((*dof).first) != 0.0)
- right_hand_side((*dof).first) = (*dof).second *
- matrix.diag_element((*dof).first);
+ new_rhs = right_hand_side((*dof).first)
+ = (*dof).second * matrix.diag_element((*dof).first);
else
{
double first_diagonal_entry = 1;
matrix.set((*dof).first, (*dof).first,
first_diagonal_entry);
- right_hand_side((*dof).first) = (*dof).second * first_diagonal_entry;
+ new_rhs = right_hand_side((*dof).first)
+ = (*dof).second * first_diagonal_entry;
};
-
+
// store the only nonzero entry
// of this line for the Gauss
// elimination step
const double diagonal_entry = matrix.diag_element((*dof).first);
-
// do the Gauss step
- for (unsigned int row=0; row<n_dofs; ++row)
+ for (unsigned int row=0; row<n_dofs; ++row)
for (int j=sparsity_rowstart[row];
- j<sparsity_rowstart[row+1]; ++row)
+ j<sparsity_rowstart[row+1]; ++j)
if ((sparsity_colnums[j] == (signed int)(*dof).first) &&
((signed int)row != (*dof).first))
// this line has an entry
{
// correct right hand side
right_hand_side(row) -= matrix.global_entry(j)/diagonal_entry *
- (*dof).second;
+ new_rhs;
// set matrix entry to zero
matrix.global_entry(j) = 0.;
void
ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
+ const FiniteElement<1> &,
map<int,double> &) const {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
Assert (false, ExcNotImplemented());
template <int dim>
void
ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
+ const FiniteElement<dim> &fe,
map<int,double> &boundary_values) const {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
-
+
+ // use two face iterators, since we need
+ // a DoF-iterator for the dof indices, but
+ // a Tria-iterator for the fe object
DoFHandler<dim>::active_face_iterator face = dof_handler->begin_active_face(),
endf = dof_handler->end_face();
+ Triangulation<dim>::active_face_iterator tface = tria->begin_active_face();
+
DirichletBC::const_iterator function_ptr;
- for (; face!=endf; ++face)
+
+ // field to store the indices of dofs
+ // initialize once to get the size right
+ // for the following fields.
+ vector<int> face_dofs;
+ face->get_dof_indices (face_dofs);
+ vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
+ vector<double> dof_values;
+
+ for (; face!=endf; ++face, ++tface)
if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
dirichlet_bc.end())
// face is subject to one of the
// get indices, physical location and
// boundary values of dofs on this
// face
- vector<int> face_dofs;
- vector<Point<dim> > dof_locations;
- vector<double> dof_values;
-
+ face_dofs.erase (face_dofs.begin(), face_dofs.end());
+ dof_values.erase (dof_values.begin(), dof_values.end());
face->get_dof_indices (face_dofs);
-
-//physical location
-// Assert (false, ExcNotImplemented());
- dof_locations.insert (dof_locations.begin(),
- face_dofs.size(), Point<dim>());
-
+ fe.face_ansatz_points (tface, dof_locations);
(*function_ptr).second->value_list (dof_locations, dof_values);
// enter into list