Write Triangulation::mesh_function.
+
Unify CellAccessor<1> and <2> by renaming
LineAccessor<dim> -> SubstructAccessor<dim,1>
QuadAccessor<dim> -> SubstructAccessor<dim,2>
and deriving CellAccessor<dim> from SubstructAccessor<dim,dim>
Do the same with DoFLineAccessor and DoFQuadAccessor
+Unify FEValues<dim> and FEFaceValues<dim> to FEValues<dim,subdim>
+ and use partial specialization.
+
+
Add support for fast assemblage of mass matrices.
Make AssemblerData a local class to Assembler again if gcc2.8 supports it.
// forward declarations
+template <int dim> class Boundary;
template <int dim> class FiniteElement;
template <int dim> class Quadrature;
* specified finite element using the
* quadrature points of the given
* quadrature rule.
+ *
+ * This function actually only fills
+ * the fields related to the unit face,
+ * the fields related to a real face (like
+ * gradients, true quadrature points, etc.)
+ * need to be initialized using the
+ * #reinit# function.
*/
FEValues (const FiniteElement<dim> &,
const Quadrature<dim> &,
* Reinitialize the gradients, Jacobi
* determinants, etc for the given cell
* and the given finite element.
+ *
+ * This function needs a boundary object
+ * passed, since this class needs to know
+ * how to handle faces which are located
+ * on the boundary of the domain. In that
+ * case, faces may be curved and the
+ * calculation of quadrature points,
+ * gradients and the like may need
+ * additional effort, depending on the
+ * mapping from the unit to the real cell
+ * (linear mappings use straight boundary
+ * segments, but higher order elements
+ * may use other ways.)
*/
void reinit (const Triangulation<dim>::cell_iterator &,
- const FiniteElement<dim> &);
+ const FiniteElement<dim> &,
+ const Boundary<dim> &);
/**
* Exception
+/**
+ Represent a finite element evaluated with a specific quadrature rule.
+ This class is an optimization which avoids evaluating the shape functions
+ at the quadrature points each time a quadrature takes place. Rather, the
+ values and gradients (and possibly higher order derivatives in future
+ versions of this library) are evaluated once and for all on the unit
+ face before doing the quadrature itself. Only the Jacobian matrix of
+ the transformation from the unit face to the real face and the integration
+ points in real space are calculated each time we move on to a new face.
+
+ The unit face is defined to be the tensor product of the interval $[0,1]$
+ in the present number of dimensions minus one. In part of the literature,
+ the convention is used that the unit cell be the tensor product of the
+ interval $[-1,1]$, which is to distinguished properly.
+
+ This class is very similar to the #FEValues# class; see there for more
+ documentation.
+ */
+template <int dim>
+class FEFaceValues {
+ public:
+
+
+
+ /**
+ * Number of quadrature points on
+ * the face.
+ */
+ const unsigned int n_quadrature_points;
+
+ /**
+ * Total number of shape functions
+ * on the cell adjacent to this face.
+ * This number is not the same as the
+ * number of shape functions of which
+ * the center is located on the face.
+ */
+ const unsigned int total_dofs;
+
+ /**
+ * Constructor. Fill all arrays with the
+ * values of the shape functions of the
+ * specified finite element using the
+ * quadrature points of the given
+ * quadrature rule for the face, which
+ * has a dimension one less than the
+ * cell.
+ *
+ * This function actually only fills
+ * the fields related to the unit face,
+ * the fields related to a real face (like
+ * gradients, true quadrature points, etc.)
+ * need to be initialized using the
+ * #reinit# function.
+ */
+ FEFaceValues (const FiniteElement<dim> &,
+ const Quadrature<dim-1> &,
+ const UpdateFields);
+
+ /**
+ * Return the value of the #i#th shape
+ * function at the #j# quadrature point.
+ */
+ double shape_value (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Return a pointer to the matrix holding
+ * all values of shape functions at all
+ * integration points, on the present cell.
+ * For the format of this matrix, see the
+ * documentation for the matrix itself.
+ */
+ const dFMatrix & get_shape_values () const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at the #j# quadrature point.
+ * If you want to get the derivative in
+ * one of the coordinate directions, use
+ * the appropriate function of the #Point#
+ * class to extract one component. Since
+ * only a reference to the gradient's value
+ * is returned, there should be no major
+ * performance drawback.
+ * The function returns the gradient on the
+ * real element, not the reference element.
+ */
+ const Point<dim> & shape_grad (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Return a pointer to the matrix holding
+ * all gradients of shape functions at all
+ * integration points, on the present cell.
+ * For the format of this matrix, see the
+ * documentation for the matrix itself.
+ */
+ const vector<vector<Point<dim> > > & get_shape_grads () const;
+
+ /**
+ * Return the position of the #i#th
+ * quadrature point in real space.
+ *
+ * For curved boundary cells, using
+ * biquadratic or higher mappings
+ * of the unit cell to the real cell,
+ * these points may not be on the
+ * plane submannifold on which the
+ * vertices of the face lie.
+ */
+ const Point<dim> & quadrature_point (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the vector of
+ * quadrature points.
+ */
+ const vector<Point<dim> > & get_quadrature_points () const;
+
+ /**
+ * Return the point in real space where
+ * the #i#th ansatz function is located
+ * (location is in the sense of where it
+ * assumes its nominal properties, e.g. at
+ * the vertex of a cell, at the center of
+ * a line, etc).
+ *
+ * This function is needed for the
+ * interpolation problem: if we want to
+ * transfer a continuous function to a
+ * finite element function by interpolation
+ * we have to take the continuous
+ * function's value at the ansatz function
+ * locations.
+ */
+ const Point<dim> & ansatz_point (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the vector of points
+ * denoting the location of the ansatz
+ * functions.
+ */
+ const vector<Point<dim> > & get_ansatz_points () const;
+
+ /**
+ * Return the Jacobi determinant times
+ * the weight of the #i#th quadrature
+ * point.
+ */
+ double JxW (const unsigned int i) const;
+
+ /**
+ * Return a pointer to the array holding
+ * the JxW values at the different
+ * quadrature points.
+ */
+ const vector<double> & get_JxW_values () const;
+
+ /**
+ * Reinitialize the gradients, Jacobi
+ * determinants, etc for the given cell
+ * and the given finite element.
+ *
+ * The constructor needs a boundary object
+ * passed, since this class needs to know
+ * how to handle faces which are located
+ * on the boundary of the domain. In that
+ * case, faces may be curved and the
+ * calculation of quadrature points,
+ * gradients and the like may need
+ * additional effort, depending on the
+ * mapping from the unit to the real cell
+ * (linear mappings use straight boundary
+ * segments, but higher order elements
+ * may use other ways.)
+ */
+ void reinit (const Triangulation<dim>::face_iterator &,
+ const FiniteElement<dim> &,
+ const Boundary<dim> &);
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidIndex,
+ int, int,
+ << "The index " << arg1
+ << " is out of range, it should be less than " << arg2);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcAccessToUninitializedField);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcCannotInitializeField);
+
+ private:
+ /**
+ * Store the values of the shape functions
+ * at the quadrature points. Rows in this
+ * matrix denote the values of a single
+ * shape function at the different points,
+ * columns are for a single point with the
+ * different shape functions.
+ */
+ dFMatrix shape_values;
+
+ /**
+ * Store the gradients of the shape
+ * functions at the quadrature points.
+ * Since unfortunately the full matrix
+ * classes of DEAL are not templated,
+ * we have to store them in an
+ * archetypic style.
+ *
+ * This field is reset each time
+ * #reinit# is called and contains the
+ * gradients on the real element, rather
+ * than on the reference element.
+ */
+ vector<vector<Point<dim> > > shape_gradients;
+
+ /**
+ * Store the gradients of the shape
+ * functions at the quadrature points on
+ * the unit cell.
+ * This field is set up upon construction
+ * of the object and contains the gradients
+ * on the reference element.
+ */
+ vector<vector<Point<dim> > > unit_shape_gradients;
+
+ /**
+ * Store an array of the weights of the
+ * quadrature points. This array is
+ * set up upon construction.
+ */
+ vector<double> weights;
+
+ /**
+ * Store an array of weights times the
+ * Jacobi determinant at the quadrature
+ * points. This function is reset each time
+ * #reinit# is called. The Jacobi determinant
+ * is actually the reciprocal value of the
+ * Jacobi matrices stored in this class,
+ * see the general documentation of this
+ * class for more information.
+ */
+ vector<double> JxW_values;
+
+ /**
+ * Array of quadrature points. This array
+ * is set up upon calling #reinit# and
+ * contains the quadrature points on the
+ * real element, rather than on the
+ * reference element.
+ */
+ vector<Point<dim> > quadrature_points;
+
+ /**
+ * Array of quadrature points in the unit
+ * cell. This array is set up upon
+ * construction and contains the quadrature
+ * points on the reference element.
+ */
+ vector<Point<dim-1> > unit_quadrature_points;
+
+ /**
+ * Array of points denoting the off-point
+ * of the ansatz functions. In real space
+ * (no-one seems to need the off-point
+ * on the unit cell, so no function is
+ * provided for this).
+ */
+ vector<Point<dim> > ansatz_points;
+
+ /**
+ * Store the jacobi matrices at the
+ * different quadrature points. This field
+ * is set each time #reinit# is called.
+ */
+ vector<dFMatrix> jacobi_matrices;
+
+ /**
+ * Store a pointer to the object describing
+ * the boundary of the domain.
+ */
+ const Boundary<dim> &boundary;
+
+ /**
+ * Store which fields are to be updated by
+ * the reinit function.
+ */
+ UpdateFields update_flags;
+};
+
+
+
+
/**
Base class for finite elements in arbitrary dimensions.
*/
vector<Point<dim> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
- const bool compute_q_points) const;
+ const bool compute_q_points,
+ const Boundary<dim> &boundary) const;
/**
* Return the ansatz points this FE has
* given face as a side. This function is
* needed for higher order elements, if
* we want to use curved boundary
- * approximations.
+ * approximations. For that reason, a
+ * boundary object has to be passed.
*
* The function assumes that the fields
* already have the right number of
* elements.
*/
virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
vector<Point<dim> > &ansatz_points) const;
/**
vector<Point<1> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
- const bool compute_q_points) const;
+ const bool compute_q_points,
+ const Boundary<1> &boundary) const;
/**
* Return the ansatz points this FE has
* At present it is not implemented.
*/
virtual void face_ansatz_points (const Triangulation<1>::face_iterator &face,
+ const Boundary<1> &boundary,
vector<Point<1> > &ansatz_points) const;
};
vector<Point<2> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<2> > &q_points,
- const bool compute_q_points) const;
+ const bool compute_q_points,
+ const Boundary<2> &boundary) const;
/**
* Return the ansatz points this FE has
* elements.
*/
virtual void face_ansatz_points (const Triangulation<2>::face_iterator &face,
+ const Boundary<2> &boundary,
vector<Point<2> > &ansatz_points) const;
};
vector<Point<dim> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
- const bool compute_q_points) const;
+ const bool compute_q_points,
+ const Boundary<dim> &boundary) const;
/**
* Return the ansatz points this FE has
* elements.
*/
virtual void face_ansatz_points (const Triangulation<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
vector<Point<dim> > &ansatz_points) const;
};
vector<Point<dim> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
- const bool compute_q_points) const;
+ const bool compute_q_points,
+ const Boundary<dim> &boundary) const;
};
vector<Point<dim> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
- const bool compute_q_points) const;
+ const bool compute_q_points,
+ const Boundary<dim> &boundary) const;
};
dVector &rhs_vector,
const Quadrature<dim> &quadrature,
const FiniteElement<dim> &fe,
- const UpdateFields &update_flags);
+ const UpdateFields &update_flags,
+ const Boundary<dim> &boundary);
/**
* Pointer to the dof handler object
* FEValues object need to be reinitialized
* on each cell.
*/
- const UpdateFields update_flags;
+ const UpdateFields update_flags;
+
+ /**
+ * Store a pointer to the object describing
+ * the boundary of the domain. This is
+ * necessary, since we may want to use
+ * curved faces of cells at the boundary
+ * when using higher order elements.
+ */
+ const Boundary<dim> &boundary;
};
* Pointer to the matrix to be assembled
* by this object.
*/
- dSMatrix &matrix;
+ dSMatrix &matrix;
/**
* Pointer to the vector to be assembled
* by this object.
*/
- dVector &rhs_vector;
+ dVector &rhs_vector;
/**
* Pointer to the finite element used for
* quadrature points.
*/
FEValues<dim> fe_values;
+
+ /**
+ * Store a pointer to the object describing
+ * the boundary of the domain. This is
+ * necessary, since we may want to use
+ * curved faces of cells at the boundary
+ * when using higher order elements.
+ */
+ const Boundary<dim> &boundary;
};
template <int dim> class Function;
template <int dim> class Equation;
template <int dim> class Assembler;
-
+template <int dim> class Boundary;
+template <int dim> class StraightBoundary;
induced by hanging nodes.
\end{itemize}
+ The #assemble# function needs an object describing the boundary of the domain,
+ since for higher order finite elements, we may be tempted to use curved faces
+ of cells for better approximation of the boundary. In this case, the
+ transformation from the unit cell to the real cell requires knowledge of
+ the exact boundary of the domain. By default it is assumed that the
+ boundary be approximated by straight segments.
+
{\bf Solving}
const Quadrature<dim> &q,
const FiniteElement<dim> &fe,
const UpdateFields &update_flags,
- const DirichletBC &dirichlet_bc = DirichletBC());
+ const DirichletBC &dirichlet_bc = DirichletBC(),
+ const Boundary<dim> &boundary = StraightBoundary<dim>());
/**
* Solve the system of equations.
dVector &difference,
const Quadrature<dim> &q,
const FiniteElement<dim> &fe,
- const NormType &norm) const;
+ const NormType &norm,
+ const Boundary<dim> &boundary=StraightBoundary<dim>()) const;
/**
* Initialize the #DataOut# object with
*/
virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc,
const FiniteElement<dim> &fe,
+ const Boundary<dim> &boundary,
map<int,double> &boundary_values) const;
/**
dVector &solution,
dVector &right_hand_side,
const DirichletBC &dirichlet_bc,
- const FiniteElement<dim> &fe);
+ const FiniteElement<dim> &fe,
+ const Boundary<dim> &boundary);
friend class Assembler<dim>;
};
#include <fe/quadrature.h>
#include <grid/tria_iterator.h>
#include <grid/tria_accessor.h>
-
+#include <grid/tria_boundary.h>
template <int dim>
double FEValues<dim>::shape_value (const unsigned int i,
const unsigned int j) const {
- Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
- Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+ Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
+ Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
return shape_values(i,j);
};
const Point<dim> &
FEValues<dim>::shape_grad (const unsigned int i,
const unsigned int j) const {
- Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
- Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+ Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
+ Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
return shape_gradients[i][j];
template <int dim>
void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
- const FiniteElement<dim> &fe) {
+ const FiniteElement<dim> &fe,
+ const Boundary<dim> &boundary) {
// fill jacobi matrices and real
// quadrature points
if ((update_flags | update_jacobians) ||
ansatz_points,
update_flags | update_ansatz_points,
quadrature_points,
- update_flags | update_q_points);
+ update_flags | update_q_points,
+ boundary);
// compute gradients on real element if
// requested
// determinant
if (update_flags | update_JxW_values)
{
- Assert (update_flags | update_jacobians,
- ExcCannotInitializeField());
+ Assert (update_flags | update_jacobians, ExcCannotInitializeField());
for (unsigned int i=0; i<n_quadrature_points; ++i)
JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
};
vector<Point<dim> > &,
const bool,
vector<Point<dim> > &,
- const bool) const {
+ const bool,
+ const Boundary<dim> &) const {
Assert (false, ExcPureFunctionCalled());
};
template <int dim>
void FiniteElementBase<dim>::face_ansatz_points (const typename Triangulation<dim>::face_iterator &,
+ const Boundary<dim> &,
vector<Point<dim> > &) const {
Assert (false, ExcPureFunctionCalled());
};
vector<Point<1> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
- const bool compute_q_points) const {
+ const bool compute_q_points,
+ const Boundary<1> &) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
Assert (q_points.size() == unit_points.size(),
void FiniteElement<1>::face_ansatz_points (const typename Triangulation<1>::face_iterator &,
+ const Boundary<1> &,
vector<Point<1> > &) const {
// is this function useful in 1D?
Assert (false, ExcPureFunctionCalled());
vector<Point<2> > &,
const bool,
vector<Point<2> > &,
- const bool) const {
+ const bool,
+ const Boundary<2> &) const {
Assert (false, ExcPureFunctionCalled());
};
void FiniteElement<2>::face_ansatz_points (const typename Triangulation<2>::face_iterator &,
+ const Boundary<2> &,
vector<Point<2> > &) const {
// is this function useful in 1D?
Assert (false, ExcPureFunctionCalled());
vector<Point<1> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
- const bool compute_q_points) const {
+ const bool compute_q_points,
+ const Boundary<1> &boundary) const {
// simply pass down
FiniteElement<1>::fill_fe_values (cell, unit_points,
jacobians, compute_jacobians,
ansatz_points, compute_ansatz_points,
- q_points, compute_q_points);
+ q_points, compute_q_points, boundary);
};
void FELinear<1>::face_ansatz_points (const Triangulation<1>::face_iterator &,
+ const Boundary<1> &,
vector<Point<1> > &) const {
Assert (false, ExcPureFunctionCalled());
};
vector<Point<2> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<2> > &q_points,
- const bool compute_q_points) const {
+ const bool compute_q_points,
+ const Boundary<2> &) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
Assert (q_points.size() == unit_points.size(),
template <int dim>
void FELinear<dim>::face_ansatz_points (const typename Triangulation<dim>::face_iterator &face,
+ const Boundary<dim> &,
vector<Point<dim> > &ansatz_points) const {
Assert (ansatz_points.size() == (1<<(dim-1)),
typename FiniteElementBase<dim>::ExcWrongFieldDimension (ansatz_points.size(),
vector<Point<1> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
- const bool compute_q_points) const {
+ const bool compute_q_points,
+ const Boundary<1> &boundary) const {
// simply pass down
FiniteElement<1>::fill_fe_values (cell, unit_points,
jacobians, compute_jacobians,
ansatz_points, compute_ansatz_points,
- q_points, compute_q_points);
+ q_points, compute_q_points, boundary);
};
vector<Point<2> > &ansatz_points,
const bool,
vector<Point<2> > &q_points,
- const bool) const {
+ const bool,
+ const Boundary<2> &) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
Assert (q_points.size() == unit_points.size(),
vector<Point<1> > &ansatz_points,
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
- const bool compute_q_points) const {
+ const bool compute_q_points,
+ const Boundary<1> &boundary) const {
// simply pass down
FiniteElement<1>::fill_fe_values (cell, unit_points,
jacobians, compute_jacobians,
ansatz_points, compute_ansatz_points,
- q_points, compute_q_points);
+ q_points, compute_q_points, boundary);
};
vector<Point<2> > &ansatz_points,
const bool,
vector<Point<2> > &q_points,
- const bool) const {
+ const bool,
+ const Boundary<2> &) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
Assert (q_points.size() == unit_points.size(),
dVector &rhs_vector,
const Quadrature<dim> &quadrature,
const FiniteElement<dim> &fe,
- const UpdateFields &update_flags) :
+ const UpdateFields &update_flags,
+ const Boundary<dim> &boundary) :
dof(dof),
assemble_matrix(assemble_matrix),
assemble_rhs(assemble_rhs),
rhs_vector(rhs_vector),
quadrature(quadrature),
fe(fe),
- update_flags(update_flags) {};
+ update_flags(update_flags),
+ boundary(boundary) {};
fe(((AssemblerData<dim>*)local_data)->fe),
fe_values (((AssemblerData<dim>*)local_data)->fe,
((AssemblerData<dim>*)local_data)->quadrature,
- ((AssemblerData<dim>*)local_data)->update_flags)
+ ((AssemblerData<dim>*)local_data)->update_flags),
+ boundary(((AssemblerData<dim>*)local_data)->boundary)
{
Assert (matrix.m() == dof_handler->n_dofs(), ExcInvalidData());
Assert (matrix.n() == dof_handler->n_dofs(), ExcInvalidData());
fe_values.reinit (Triangulation<dim>::cell_iterator (tria,
present_level,
present_index),
- fe);
+ fe,
+ boundary);
const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs;
// clear cell matrix
const Quadrature<dim> &quadrature,
const FiniteElement<dim> &fe,
const UpdateFields &update_flags,
- const DirichletBC &dirichlet_bc) {
+ const DirichletBC &dirichlet_bc,
+ const Boundary<dim> &boundary) {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
system_sparsity.reinit (dof_handler->n_dofs(),
right_hand_side,
quadrature,
fe,
- update_flags);
+ update_flags,
+ boundary);
active_assemble_iterator assembler (tria,
tria->begin_active()->level(),
tria->begin_active()->index(),
// in the docs
apply_dirichlet_bc (system_matrix, solution,
right_hand_side,
- dirichlet_bc, fe);
+ dirichlet_bc, fe, boundary);
};
dVector &difference,
const Quadrature<dim> &q,
const FiniteElement<dim> &fe,
- const NormType &norm) const {
+ const NormType &norm,
+ const Boundary<dim> &boundary) const {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE());
{
double diff=0;
// initialize for this cell
- fe_values.reinit (tria_cell, fe);
+ fe_values.reinit (tria_cell, fe, boundary);
switch (norm)
{
dVector &solution,
dVector &right_hand_side,
const DirichletBC &dirichlet_bc,
- const FiniteElement<dim> &fe) {
+ const FiniteElement<dim> &fe,
+ const Boundary<dim> &boundary) {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
ExcInvalidBoundaryIndicator());
// 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, fe, boundary_values);
+ make_boundary_value_list (dirichlet_bc, fe, boundary, boundary_values);
map<int,double>::const_iterator dof, endd;
const unsigned int n_dofs = (unsigned int)matrix.m();
void
ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
const FiniteElement<1> &,
+ const Boundary<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,
+ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
const FiniteElement<dim> &fe,
+ const Boundary<dim> &boundary,
map<int,double> &boundary_values) const {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
face_dofs.erase (face_dofs.begin(), face_dofs.end());
dof_values.erase (dof_values.begin(), dof_values.end());
face->get_dof_indices (face_dofs);
- fe.face_ansatz_points (tface, dof_locations);
+ fe.face_ansatz_points (tface, boundary, dof_locations);
(*function_ptr).second->value_list (dof_locations, dof_values);
// enter into list