/*---------------------------- function.h ---------------------------*/
-//forward declaration
-template <int dim> class Point;
-template <class T, class Alloc=alloc> class vector;
+#include <base/exceptions.h>
+#include <vector.h>
+#include <grid/point.h>
+
+/**
+ This class is a model for a continuous function. It returns the value
+ of the function at a given point through the #operator ()# member function,
+ which is virtual. It also has a function to return a whole list of function
+ values at different points to reduce the overhead of the virtual function
+ calls; this function is preset to successively call the function returning
+ one value at a time.
+
+ There are other functions return the gradient of the function at one or
+ several points. You only have to overload those functions you need; the
+ functions returning several values at a time will call those returning
+ only one value, while those ones will throw an exception when called but
+ not overloaded.
+
+ Unless only called a very small number of times, you should overload
+ both those functions returning only one value as well as those returning
+ a whole array, since the cost of evaluation of a point value is often
+ less than the virtual function call itself.
+*/
template <int dim>
class Function {
public:
- virtual double operator () (const Point<dim> &p) const = 0;
+ /**
+ * Return the value of the function
+ * at the given point.
+ */
+ virtual double operator () (const Point<dim> &p) const;
+
+ /**
+ * Set #values# to the point values
+ * of the function at the #points#.
+ * It is assumed that #values# be
+ * empty.
+ */
virtual void value_list (const vector<Point<dim> > &points,
vector<double> &values) const;
+
+ /**
+ * Return the gradient of the function
+ * at the given point.
+ */
+ virtual Point<dim> gradient (const Point<dim> &p) const;
+
+ /**
+ * Set #gradients# to the gradients of
+ * the function at the #points#.
+ * It is assumed that #values# be
+ * empty.
+ */
+ virtual void gradient_list (const vector<Point<dim> > &points,
+ vector<Point<dim> > &gradients) const;
+
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcPureFunctionCalled);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcVectorNotEmpty);
};
#include <grid/tria_boundary.h>
#include <grid/dof_constraints.h>
#include <basic/data_io.h>
+#include <basic/function.h>
#include <fe/fe_lib.h>
#include <fe/quadrature_lib.h>
#include <numerics/base.h>
#include <numerics/assembler.h>
#include <lac/dsmatrix.h>
+
+#include <map.h>
#include <fstream.h>
#include <cmath>
extern "C" {
<< "Invalid index " << arg1
<< ", index must be between " << arg2
<< " and " << arg3 << ".");
+ /**
+ * Exception
+ */
+ DeclException0 (ExcVectorNotEmpty);
protected:
/**
const unsigned int i,
const int index) const;
+ /**
+ * Return the indices of the dofs of this
+ * line in the standard ordering: dofs
+ * on vertex 0, dofs on vertex 1,
+ * dofs on line.
+ *
+ * The dof indices are appended to the end
+ * of the vector. It is assumed that
+ * the vector be empty beforehand.
+ */
+ void get_dof_indices (vector<int> &dof_indices) const;
+
/**
* Return the #i#th child as a DoF line
* iterator. This function is needed since
const unsigned int i,
const int index) const;
+ /**
+ * Return the indices of the dofs of this
+ * quad in the standard ordering: dofs
+ * on vertex 0, dofs on vertex 1, etc,
+ * dofs on line 0, dofs on line 1, etc,
+ * dofs on quad 0, etc.
+ *
+ * The dof indices are appended to the end
+ * of the vector. It is assumed that
+ * the vector be empty beforehand.
+ */
+ void get_dof_indices (vector<int> &dof_indices) const;
+
/**
* Return a pointer to the #i#th line
* bounding this #Quad#.
DoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {};
// do this here, since this way the
// CellAccessor has the possibility to know
- // what a substruct_iterator is. Otherwise
+ // what a face_iterator is. Otherwise
// it would have to ask the DoFHandler<dim>
// which would need another big include
// file and maybe cyclic interdependence
- typedef void * substruct_iterator;
+ typedef void * face_iterator;
};
DoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {};
// do this here, since this way the
// CellAccessor has the possibility to know
- // what a substruct_iterator is. Otherwise
+ // what a face_iterator is. Otherwise
// it would have to ask the DoFHandler<dim>
// which would need another big include
// file and maybe cyclic interdependence
- typedef TriaIterator<2,DoFLineAccessor<2,LineAccessor<2> > > substruct_iterator;
+ typedef TriaIterator<2,DoFLineAccessor<2,LineAccessor<2> > > face_iterator;
};
* This function is not implemented in 1D,
* and maps to DoFQuadAccessor::line in 2D.
*/
- typename DoFSubstructAccessor<dim>::substruct_iterator
+ typename DoFSubstructAccessor<dim>::face_iterator
face (const unsigned int i) const;
- /**
- * Return the indices of the dofs of this
- * cell in the standard ordering: dofs
- * on vertex 0, dofs on vertex 1, etc,
- * dofs on line 0, dofs on line 1, etc,
- * dofs on quad 0, etc.
- *
- * The dof indices are appended to the end
- * of the vector. It is assumed that
- * the vector be empty beforehand.
- */
- void get_dof_indices (vector<int> &dof_indices) const;
-
/**
* Return the value of the given vector
* restricted to the dofs of this
void get_dof_values (const dVector &values,
vector<double> &dof_values) const;
- /**
- * Exception
- */
- DeclException0 (ExcVectorNotEmpty);
/**
* Exception
*/
typedef line_iterator cell_iterator;
typedef active_line_iterator active_cell_iterator;
- typedef void * raw_substruct_iterator;
- typedef void * substruct_iterator;
- typedef void * active_substruct_iterator;
+ typedef void * raw_face_iterator;
+ typedef void * face_iterator;
+ typedef void * active_face_iterator;
};
typedef quad_iterator cell_iterator;
typedef active_quad_iterator active_cell_iterator;
- typedef raw_line_iterator raw_substruct_iterator;
- typedef line_iterator substruct_iterator;
- typedef active_line_iterator active_substruct_iterator;
+ typedef raw_line_iterator raw_face_iterator;
+ typedef line_iterator face_iterator;
+ typedef active_line_iterator active_face_iterator;
};
typedef typename DoFDimensionInfo<dim>::cell_iterator cell_iterator;
typedef typename DoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
- typedef typename DoFDimensionInfo<dim>::raw_substruct_iterator raw_substruct_iterator;
- typedef typename DoFDimensionInfo<dim>::substruct_iterator substruct_iterator;
- typedef typename DoFDimensionInfo<dim>::active_substruct_iterator active_substruct_iterator;
+ typedef typename DoFDimensionInfo<dim>::raw_face_iterator raw_face_iterator;
+ typedef typename DoFDimensionInfo<dim>::face_iterator face_iterator;
+ typedef typename DoFDimensionInfo<dim>::active_face_iterator active_face_iterator;
/**
/*---------------------------------------*/
+ /**
+ * @name Face iterator functions
+ */
+ /*@{*/
+ /**
+ * Return iterator to the first face, used
+ * or not, on level #level#. If a level
+ * has no faces, a past-the-end iterator
+ * is returned.
+ *
+ * This function calls #begin_raw_line#
+ * in 2D and #begin_raw_quad# in 3D.
+ */
+ raw_face_iterator begin_raw_face (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first used face
+ * on level #level#.
+ *
+ * This function calls #begin_line#
+ * in 2D and #begin_quad# in 3D.
+ */
+ face_iterator begin_face (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first active
+ * face on level #level#.
+ *
+ * This function calls #begin_active_line#
+ * in 2D and #begin_active_quad# in 3D.
+ */
+ active_face_iterator begin_active_face(const unsigned int level = 0) const;
+
+ /**
+ * Return iterator past the end; this
+ * iterator serves for comparisons of
+ * iterators with past-the-end or
+ * before-the-beginning states.
+ *
+ * This function calls #end_line#
+ * in 2D and #end_quad# in 3D.
+ */
+ raw_face_iterator end_face () const;
+
+ /**
+ * Return an iterator pointing to the
+ * last face, used or not.
+ *
+ * This function calls #last_raw_line#
+ * in 2D and #last_raw_quad# in 3D.
+ */
+ raw_face_iterator last_raw_face () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * face of the level #level#, used or not.
+ *
+ * This function calls #last_raw_line#
+ * in 2D and #last_raw_quad# in 3D.
+ */
+ raw_face_iterator last_raw_face (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used face.
+ *
+ * This function calls #last_line#
+ * in 2D and #last_quad# in 3D.
+ */
+ face_iterator last_face () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used face on level #level#.
+ *
+ * This function calls #last_line#
+ * in 2D and #last_quad# in 3D.
+ */
+ face_iterator last_face (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active face.
+ *
+ * This function calls #last_active_line#
+ * in 2D and #last_active_quad# in 3D.
+ */
+ active_face_iterator last_active_face () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active face on level #level#.
+ *
+ * This function calls #last_active_line#
+ * in 2D and #last_active_quad# in 3D.
+ */
+ active_face_iterator last_active_face (const unsigned int level) const;
+ //@}
+
+
+ /*---------------------------------------*/
+
/**
* @name Line iterator functions
*/
DeclException1 (ExcMatrixHasWrongSize,
int,
<< "The matrix has the wrong dimension " << arg1);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcFunctionNotUseful);
protected:
/**
#void *#. Thus these types exist, but are useless and will
certainly make any involuntary use visible.
- The same applies for the #substruct_iterator# types, since lines
+ The same applies for the #face_iterator# types, since lines
have no substructures apart from vertices, which are handled in
a different way, however.
*/
typedef line_iterator cell_iterator;
typedef active_line_iterator active_cell_iterator;
- typedef void * raw_substruct_iterator;
- typedef void * substruct_iterator;
- typedef void * active_substruct_iterator;
+ typedef void * raw_face_iterator;
+ typedef void * face_iterator;
+ typedef void * active_face_iterator;
};
typedef quad_iterator cell_iterator;
typedef active_quad_iterator active_cell_iterator;
- typedef raw_line_iterator raw_substruct_iterator;
- typedef line_iterator substruct_iterator;
- typedef active_line_iterator active_substruct_iterator;
+ typedef raw_line_iterator raw_face_iterator;
+ typedef line_iterator face_iterator;
+ typedef active_line_iterator active_face_iterator;
\end{verbatim}
*/
class TriaDimensionInfo<2> {
typedef quad_iterator cell_iterator;
typedef active_quad_iterator active_cell_iterator;
- typedef raw_line_iterator raw_substruct_iterator;
- typedef line_iterator substruct_iterator;
- typedef active_line_iterator active_substruct_iterator;
+ typedef raw_line_iterator raw_face_iterator;
+ typedef line_iterator face_iterator;
+ typedef active_line_iterator active_face_iterator;
};
typedef quad_iterator cell_iterator;
typedef active_quad_iterator active_cell_iterator;
- typedef raw_line_iterator raw_substruct_iterator;
- typedef line_iterator substruct_iterator;
- typedef active_line_iterator active_substruct_iterator;
+ typedef raw_line_iterator raw_face_iterator;
+ typedef line_iterator face_iterator;
+ typedef active_line_iterator active_face_iterator;
\end{verbatim}
By using the cell iterators, you can write code nearly independent of
typedef typename TriaDimensionInfo<dim>::cell_iterator cell_iterator;
typedef typename TriaDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
- typedef typename TriaDimensionInfo<dim>::raw_substruct_iterator raw_substruct_iterator;
- typedef typename TriaDimensionInfo<dim>::substruct_iterator substruct_iterator;
- typedef typename TriaDimensionInfo<dim>::active_substruct_iterator active_substruct_iterator;
+ typedef typename TriaDimensionInfo<dim>::raw_face_iterator raw_face_iterator;
+ typedef typename TriaDimensionInfo<dim>::face_iterator face_iterator;
+ typedef typename TriaDimensionInfo<dim>::active_face_iterator active_face_iterator;
/**
* Create a triangulation and create
//@}
/*---------------------------------------*/
+ /*---------------------------------------*/
+
+ /**
+ * @name Face iterator functions
+ */
+ /*@{*/
+ /**
+ * Return iterator to the first face, used
+ * or not, on level #level#. If a level
+ * has no faces, a past-the-end iterator
+ * is returned.
+ *
+ * This function calls #begin_raw_line#
+ * in 2D and #begin_raw_quad# in 3D.
+ */
+ raw_face_iterator begin_raw_face (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first used face
+ * on level #level#.
+ *
+ * This function calls #begin_line#
+ * in 2D and #begin_quad# in 3D.
+ */
+ face_iterator begin_face (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first active
+ * face on level #level#.
+ *
+ * This function calls #begin_active_line#
+ * in 2D and #begin_active_quad# in 3D.
+ */
+ active_face_iterator begin_active_face(const unsigned int level = 0) const;
+
+ /**
+ * Return iterator past the end; this
+ * iterator serves for comparisons of
+ * iterators with past-the-end or
+ * before-the-beginning states.
+ *
+ * This function calls #end_line#
+ * in 2D and #end_quad# in 3D.
+ */
+ raw_face_iterator end_face () const;
+
+ /**
+ * Return an iterator pointing to the
+ * last face, used or not.
+ *
+ * This function calls #last_raw_line#
+ * in 2D and #last_raw_quad# in 3D.
+ */
+ raw_face_iterator last_raw_face () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * face of the level #level#, used or not.
+ *
+ * This function calls #last_raw_line#
+ * in 2D and #last_raw_quad# in 3D.
+ */
+ raw_face_iterator last_raw_face (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used face.
+ *
+ * This function calls #last_line#
+ * in 2D and #last_quad# in 3D.
+ */
+ face_iterator last_face () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used face on level #level#.
+ *
+ * This function calls #last_line#
+ * in 2D and #last_quad# in 3D.
+ */
+ face_iterator last_face (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active face.
+ *
+ * This function calls #last_active_line#
+ * in 2D and #last_active_quad# in 3D.
+ */
+ active_face_iterator last_active_face () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active face on level #level#.
+ *
+ * This function calls #last_active_line#
+ * in 2D and #last_active_quad# in 3D.
+ */
+ active_face_iterator last_active_face (const unsigned int level) const;
+ //@}
+
+
+ /*---------------------------------------*/
/**
* @name Line iterator functions
template <int dim> class Quadrature;
template <int dim> class DataOut;
template <int dim> class Function;
+template <int dim> class Equation;
+template <int dim> class Assembler;
template <class Key, class T, class Compare = less<Key>, class Alloc = alloc> class map;
eliminating all coupling between this degree of freedom and others. Now
also the column consists only of zeroes, apart from the main diagonal entry.
+ It seems as if we had to make clear not to overwrite the lines of other
+ boundary nodes when doing the Gauss elimination step. However, since we
+ reset the right hand side when passing such a node, it is not a problem
+ to change the right hand side values of other boundary nodes not yet
+ processed. It would be a problem to change those entries of nodes already
+ processed, but since the matrix entry of the present column on the row
+ of an already processed node is zero, the Gauss step does not change
+ the right hand side. We need therefore not take special care of other
+ boundary nodes.
+
To make solving faster, we preset the solution vector with the right boundary
values. Since boundary nodes can never be hanging nodes, and since all other
entries of the solution vector are zero, we need not condense the solution
since their search directions may contains components in the direction
of the boundary node. For this reason, we perform a very simple line
balancing by not setting the main diagonal entry to unity, but rather
- to the value it had before deleting this line. Of course we have to change
+ to the value it had before deleting this line, or to the first nonzero
+ main diagonal entry if it is zero from a previous Gauss elimination
+ step. Of course we have to change
the right hand side appropriately. This is not a very good
strategy, but it at least should give the main diagonal entry a value
in the right order of dimension, which makes the solving process a bit
not during solving iteratively, it may or may not be necessary to set
the correct value after solving again. This question is an open one as of
now and may be answered by future experience.
+
+ At present, boundary values are interpolated, i.e. a node is given the
+ point value of the boundary function. In some cases, it may be necessary
+ to use the L2-projection of the boundary function or any other method.
+ This can be done by overloading the virtual function
+ #make_boundary_value_list# which must return a list of boundary dofs
+ and their corresponding values.
{\bf Computing errors}
* See the general documentation of this
* class for more detail.
*/
- typedef map<unsigned char,Function<dim> > DirichletBC;
+ typedef map<unsigned char,Function<dim>* > DirichletBC;
/**
* Constructor. Use this triangulation and
* want anything else.
*/
virtual pair<char*,char*> get_solution_name () const;
+
+ /**
+ * Make up the list of node subject
+ * to Dirichlet boundary conditions
+ * and the values they are to be
+ * assigned.
+ *
+ * See the general doc for more
+ * information.
+ */
+ virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc,
+ map<int,double> &boundary_values) const;
/**
* Exception
* Exception
*/
DeclException0 (ExcNotImplemented);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidBoundaryIndicator);
protected:
/**
* as described in the general
* documentation.
*/
- void apply_dirichlet_bc (dSMatrix &matrix,
- dVector &solution,
- dVector &right_hand_side,
+ void apply_dirichlet_bc (dSMatrix &matrix,
+ dVector &solution,
+ dVector &right_hand_side,
const DirichletBC &dirichlet_bc);
friend class Assembler<dim>;
+template <int dim, class BaseClass>
+void
+DoFLineAccessor<dim,BaseClass>::get_dof_indices (vector<int> &dof_indices) const {
+ Assert (dof_handler != 0, ExcInvalidObject());
+ Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+ Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
+
+ dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
+ for (unsigned int vertex=0; vertex<2; ++vertex)
+ for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
+ dof_indices.push_back (vertex_dof_index(vertex,d));
+ for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
+ dof_indices.push_back (dof_index(d));
+};
+
+
+
template <int dim, class BaseClass>
TriaIterator<dim,DoFLineAccessor<dim,BaseClass> >
DoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
+template <int dim, class BaseClass>
+void
+DoFQuadAccessor<dim,BaseClass>::get_dof_indices (vector<int> &dof_indices) const {
+ Assert (dof_handler != 0, ExcInvalidObject());
+ Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+ Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
+
+ dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
+ dof_indices.push_back (vertex_dof_index(vertex,d));
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
+ dof_indices.push_back (this->line(line)->dof_index(d));
+ for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_quad; ++d)
+ dof_indices.push_back (dof_index(d));
+};
+
+
+
+
template <int dim, class BaseClass>
TriaIterator<dim,DoFLineAccessor<dim,LineAccessor<dim> > >
DoFQuadAccessor<dim,BaseClass>::line (const unsigned int i) const {
-DoFSubstructAccessor<1>::substruct_iterator
+DoFSubstructAccessor<1>::face_iterator
DoFCellAccessor<1>::face (const unsigned int) const {
Assert (false, ExcNotUsefulForThisDimension());
return 0;
-DoFSubstructAccessor<2>::substruct_iterator
+DoFSubstructAccessor<2>::face_iterator
DoFCellAccessor<2>::face (const unsigned int i) const {
return line(i);
};
-void
-DoFCellAccessor<1>::get_dof_indices (vector<int> &dof_indices) const {
- Assert (dof_handler != 0, ExcInvalidObject());
- Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
- Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
-
- dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
- for (unsigned int vertex=0; vertex<2; ++vertex)
- for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
- dof_indices.push_back (vertex_dof_index(vertex,d));
- for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
- dof_indices.push_back (dof_index(d));
-};
-
-
-
-
-void
-DoFCellAccessor<2>::get_dof_indices (vector<int> &dof_indices) const {
- Assert (dof_handler != 0, ExcInvalidObject());
- Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
- Assert (dof_indices.size() == 0, ExcVectorNotEmpty());
-
- dof_indices.reserve (dof_handler->get_selected_fe().total_dofs);
- for (unsigned int vertex=0; vertex<4; ++vertex)
- for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_vertex; ++d)
- dof_indices.push_back (vertex_dof_index(vertex,d));
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_line; ++d)
- dof_indices.push_back (this->line(line)->dof_index(d));
- for (unsigned int d=0; d<dof_handler->get_selected_fe().dofs_per_quad; ++d)
- dof_indices.push_back (dof_index(d));
-};
-
-
-
-
void
DoFCellAccessor<1>::get_dof_values (const dVector &values,
vector<double> &dof_values) const {
void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed,
dSMatrixStruct &condensed) const {
Assert (sorted == true, ExcMatrixNotClosed());
- Assert (uncondensed.compressed == true, ExcMatrixNotClosed());
+ Assert (uncondensed.is_compressed() == true, ExcMatrixNotClosed());
Assert (uncondensed.n_rows() == uncondensed.n_cols(),
ExcMatrixNotSquare());
next_constraint = lines.begin();
- for (int row=0; row<uncondensed.rows; ++row)
+ for (int row=0; row<uncondensed.n_rows(); ++row)
if (new_line[row] != -1)
// line not constrained
// copy entries if column will not
// be condensed away, distribute
// otherwise
- for (int j=uncondensed.rowstart[row]; j<uncondensed.rowstart[row+1]; ++j)
- if (new_line[uncondensed.colnums[j]] != -1)
- condensed.add (new_line[row], new_line[uncondensed.colnums[j]]);
+ for (int j=uncondensed.get_rowstart_indices()[row];
+ j<uncondensed.get_rowstart_indices()[row+1]; ++j)
+ if (new_line[uncondensed.get_column_numbers()[j]] != -1)
+ condensed.add (new_line[row], new_line[uncondensed.get_column_numbers()[j]]);
else
{
// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while ((*c).line != (unsigned int)uncondensed.colnums[j]) ++c;
+ while ((*c).line != (unsigned int)uncondensed.get_column_numbers()[j])
+ ++c;
for (unsigned int q=0; q!=(*c).entries.size(); ++q)
condensed.add (new_line[row], new_line[(*c).entries[q].first]);
else
// line must be distributed
{
- for (int j=uncondensed.rowstart[row]; j<uncondensed.rowstart[row+1]; ++j)
+ for (int j=uncondensed.get_rowstart_indices()[row];
+ j<uncondensed.get_rowstart_indices()[row+1]; ++j)
// for each entry: distribute
- if (new_line[uncondensed.colnums[j]] != -1)
+ if (new_line[uncondensed.get_column_numbers()[j]] != -1)
// column is not constrained
for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
condensed.add (new_line[(*next_constraint).entries[q].first],
- new_line[uncondensed.colnums[j]]);
+ new_line[uncondensed.get_column_numbers()[j]]);
else
// not only this line but
// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while ((*c).line != (unsigned int)uncondensed.colnums[j]) ++c;
+ while ((*c).line != (unsigned int)uncondensed.get_column_numbers()[j]) ++c;
for (unsigned int p=0; p!=(*c).entries.size(); ++p)
for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const {
Assert (sorted == true, ExcMatrixNotClosed());
- Assert (sparsity.compressed == false, ExcMatrixIsClosed());
+ Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
Assert (sparsity.n_rows() == sparsity.n_cols(),
ExcMatrixNotSquare());
for (int row=0; row<n_rows; ++row)
if (distribute[row] == -1)
// regular line. loop over cols
- for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+ for (int j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
// end of row reached?
- if (sparsity.colnums[j] == -1)
+ if (sparsity.get_column_numbers()[j] == -1)
break;
else
{
- if (distribute[sparsity.colnums[j]] != -1)
+ if (distribute[sparsity.get_column_numbers()[j]] != -1)
// distribute entry at regular
// row #row# and irregular column
// sparsity.colnums[j]
for (unsigned int q=0;
- q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+ q!=lines[distribute[sparsity.get_column_numbers()[j]]].entries.size();
+ ++q)
sparsity.add (row,
- lines[distribute[sparsity.colnums[j]]].entries[q].first);
+ lines[distribute[sparsity.get_column_numbers()[j]]].
+ entries[q].first);
}
else
// row must be distributed
- for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+ for (int j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
// end of row reached?
- if (sparsity.colnums[j] == -1)
+ if (sparsity.get_column_numbers()[j] == -1)
break;
else
{
- if (distribute[sparsity.colnums[j]] == -1)
+ if (distribute[sparsity.get_column_numbers()[j]] == -1)
// distribute entry at irregular
// row #row# and regular column
// sparsity.colnums[j]
for (unsigned int q=0;
q!=lines[distribute[row]].entries.size(); ++q)
sparsity.add (lines[distribute[row]].entries[q].first,
- sparsity.colnums[j]);
+ sparsity.get_column_numbers()[j]);
else
// distribute entry at irregular
// row #row# and irregular column
- // sparsity.colnums[j]
+ // sparsity.get_column_numbers()[j]
for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
for (unsigned int q=0;
- q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+ q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries.size(); ++q)
sparsity.add (lines[distribute[row]].entries[p].first,
- lines[distribute[sparsity.colnums[j]]].entries[q].first);
+ lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries[q].first);
};
sparsity.compress();
};
const dSMatrixStruct &uncondensed_struct = uncondensed.get_sparsity_pattern ();
Assert (sorted == true, ExcMatrixNotClosed());
- Assert (uncondensed_struct.compressed == true, ExcMatrixNotClosed());
- Assert (condensed.get_sparsity_pattern().compressed == true, ExcMatrixNotClosed());
+ Assert (uncondensed_struct.is_compressed() == true, ExcMatrixNotClosed());
+ Assert (condensed.get_sparsity_pattern().is_compressed() == true, ExcMatrixNotClosed());
Assert (uncondensed_struct.n_rows() == uncondensed_struct.n_cols(),
ExcMatrixNotSquare());
Assert (condensed.n() == condensed.m(),
next_constraint = lines.begin();
- for (int row=0; row<uncondensed_struct.rows; ++row)
+ for (int row=0; row<uncondensed_struct.n_rows(); ++row)
if (new_line[row] != -1)
// line not constrained
// copy entries if column will not
// be condensed away, distribute
// otherwise
- for (int j=uncondensed_struct.rowstart[row]; j<uncondensed_struct.rowstart[row+1]; ++j)
- if (new_line[uncondensed_struct.colnums[j]] != -1)
- condensed.add (new_line[row], new_line[uncondensed_struct.colnums[j]],
- uncondensed.val[j]);
+ for (int j=uncondensed_struct.get_rowstart_indices()[row];
+ j<uncondensed_struct.get_rowstart_indices()[row+1]; ++j)
+ if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
+ condensed.add (new_line[row], new_line[uncondensed_struct.get_column_numbers()[j]],
+ uncondensed.global_entry(j));
else
{
// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c;
+ while ((*c).line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) ++c;
for (unsigned int q=0; q!=(*c).entries.size(); ++q)
// distribute to rows with
// appropriate weight
condensed.add (new_line[row], new_line[(*c).entries[q].first],
- uncondensed.val[j] * (*c).entries[q].second);
+ uncondensed.global_entry(j) * (*c).entries[q].second);
}
else
// line must be distributed
{
- for (int j=uncondensed_struct.rowstart[row]; j<uncondensed_struct.rowstart[row+1]; ++j)
+ for (int j=uncondensed_struct.get_rowstart_indices()[row];
+ j<uncondensed_struct.get_rowstart_indices()[row+1]; ++j)
// for each column: distribute
- if (new_line[uncondensed_struct.colnums[j]] != -1)
+ if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
// column is not constrained
for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
condensed.add (new_line[(*next_constraint).entries[q].first],
- new_line[uncondensed_struct.colnums[j]],
- uncondensed.val[j] * (*next_constraint).entries[q].second);
+ new_line[uncondensed_struct.get_column_numbers()[j]],
+ uncondensed.global_entry(j) *
+ (*next_constraint).entries[q].second);
else
// not only this line but
// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c;
+ while ((*c).line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) ++c;
for (unsigned int p=0; p!=(*c).entries.size(); ++p)
for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
condensed.add (new_line[(*next_constraint).entries[q].first],
new_line[(*c).entries[p].first],
- uncondensed.val[j] *
+ uncondensed.global_entry(j) *
(*next_constraint).entries[q].second *
(*c).entries[p].second);
};
const dSMatrixStruct &sparsity = uncondensed.get_sparsity_pattern ();
Assert (sorted == true, ExcMatrixNotClosed());
- Assert (sparsity.compressed == true, ExcMatrixNotClosed());
+ Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed());
Assert (sparsity.n_rows() == sparsity.n_cols(),
ExcMatrixNotSquare());
if (distribute[row] == -1)
// regular line. loop over cols
- for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+ for (int j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
// end of row reached?
- if (sparsity.colnums[j] == -1)
+ if (sparsity.get_column_numbers()[j] == -1)
break;
else
{
- if (distribute[sparsity.colnums[j]] != -1)
+ if (distribute[sparsity.get_column_numbers()[j]] != -1)
// distribute entry at regular
// row #row# and irregular column
- // sparsity.colnums[j]; set old
+ // sparsity.get_column_numbers()[j]; set old
// entry to zero
{
for (unsigned int q=0;
- q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+ q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries.size(); ++q)
uncondensed.add (row,
- lines[distribute[sparsity.colnums[j]]].entries[q].first,
- uncondensed.val[j] *
- lines[distribute[sparsity.colnums[j]]].entries[q].second);
- uncondensed.val[j] = 0.;
+ lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries[q].first,
+ uncondensed.global_entry(j) *
+ lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries[q].second);
+ uncondensed.global_entry(j) = 0.;
};
}
else
// row must be distributed
- for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+ for (int j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
// end of row reached?
- if (sparsity.colnums[j] == -1)
+ if (sparsity.get_column_numbers()[j] == -1)
break;
else
{
- if (distribute[sparsity.colnums[j]] == -1)
+ if (distribute[sparsity.get_column_numbers()[j]] == -1)
// distribute entry at irregular
// row #row# and regular column
- // sparsity.colnums[j]. set old
+ // sparsity.get_column_numbers()[j]. set old
// entry to zero
{
for (unsigned int q=0;
q!=lines[distribute[row]].entries.size(); ++q)
uncondensed.add (lines[distribute[row]].entries[q].first,
- sparsity.colnums[j],
- uncondensed.val[j] *
+ sparsity.get_column_numbers()[j],
+ uncondensed.global_entry(j) *
lines[distribute[row]].entries[q].second);
- uncondensed.val[j] = 0.;
+ uncondensed.global_entry(j) = 0.;
}
else
// distribute entry at irregular
// row #row# and irregular column
- // sparsity.colnums[j]
+ // sparsity.get_column_numbers()[j]
// set old entry to one if on main
// diagonal, zero otherwise
{
for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
for (unsigned int q=0;
- q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+ q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries.size(); ++q)
uncondensed.add (lines[distribute[row]].entries[p].first,
- lines[distribute[sparsity.colnums[j]]].entries[q].first,
- uncondensed.val[j] *
+ lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries[q].first,
+ uncondensed.global_entry(j) *
lines[distribute[row]].entries[p].second *
- lines[distribute[sparsity.colnums[j]]].entries[q].second);
+ lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries[q].second);
- uncondensed.val[j] = (row == sparsity.colnums[j] ?
+ uncondensed.global_entry(j) = (row == sparsity.get_column_numbers()[j] ?
1. : 0. );
};
};
//------------------------------------------------------------------
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::begin_raw_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::face_iterator
+DoFHandler<1>::begin_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::active_face_iterator
+DoFHandler<1>::begin_active_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::end_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::last_raw_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::raw_face_iterator
+DoFHandler<1>::last_raw_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::face_iterator
+DoFHandler<1>::last_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::face_iterator
+DoFHandler<1>::last_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::active_face_iterator
+DoFHandler<1>::last_active_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+DoFDimensionInfo<1>::active_face_iterator
+DoFHandler<1>::last_active_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::begin_raw_face (const unsigned int level) const {
+ return begin_raw_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::face_iterator
+DoFHandler<2>::begin_face (const unsigned int level) const {
+ return begin_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::active_face_iterator
+DoFHandler<2>::begin_active_face (const unsigned int level) const {
+ return begin_active_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::end_face () const {
+ return end_line ();
+};
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::last_raw_face () const {
+ return last_raw_line ();
+};
+
+
+
+DoFDimensionInfo<2>::raw_face_iterator
+DoFHandler<2>::last_raw_face (const unsigned int level) const {
+ return last_raw_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::face_iterator
+DoFHandler<2>::last_face () const {
+ return last_line ();
+};
+
+
+
+DoFDimensionInfo<2>::face_iterator
+DoFHandler<2>::last_face (const unsigned int level) const {
+ return last_line (level);
+};
+
+
+
+DoFDimensionInfo<2>::active_face_iterator
+DoFHandler<2>::last_active_face () const {
+ return last_active_line ();
+};
+
+
+
+DoFDimensionInfo<2>::active_face_iterator
+DoFHandler<2>::last_active_face (const unsigned int level) const {
+ return last_active_line (level);
+};
+
+
+
+// ---------------------------------------------------------------
+
+
+
template <int dim>
typename DoFHandler<dim>::raw_line_iterator
DoFHandler<dim>::begin_raw_line (const unsigned int level) const {
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::begin_raw_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::face_iterator
+Triangulation<1>::begin_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::active_face_iterator
+Triangulation<1>::begin_active_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::end_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::last_raw_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::raw_face_iterator
+Triangulation<1>::last_raw_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::face_iterator
+Triangulation<1>::last_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::face_iterator
+Triangulation<1>::last_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::active_face_iterator
+Triangulation<1>::last_active_face () const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+TriaDimensionInfo<1>::active_face_iterator
+Triangulation<1>::last_active_face (const unsigned int) const {
+ Assert (false, ExcFunctionNotUseful());
+ return 0;
+};
+
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::begin_raw_face (const unsigned int level) const {
+ return begin_raw_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::face_iterator
+Triangulation<2>::begin_face (const unsigned int level) const {
+ return begin_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::active_face_iterator
+Triangulation<2>::begin_active_face (const unsigned int level) const {
+ return begin_active_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::end_face () const {
+ return end_line ();
+};
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::last_raw_face () const {
+ return last_raw_line ();
+};
+
+
+
+TriaDimensionInfo<2>::raw_face_iterator
+Triangulation<2>::last_raw_face (const unsigned int level) const {
+ return last_raw_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::face_iterator
+Triangulation<2>::last_face () const {
+ return last_line ();
+};
+
+
+
+TriaDimensionInfo<2>::face_iterator
+Triangulation<2>::last_face (const unsigned int level) const {
+ return last_line (level);
+};
+
+
+
+TriaDimensionInfo<2>::active_face_iterator
+Triangulation<2>::last_active_face () const {
+ return last_active_line ();
+};
+
+
+
+TriaDimensionInfo<2>::active_face_iterator
+Triangulation<2>::last_active_face (const unsigned int level) const {
+ return last_active_line (level);
+};
+
+
+
+
+
template <int dim>
typename TriaDimensionInfo<dim>::raw_line_iterator
Triangulation<dim>::begin_raw_line (unsigned int level) const {
-Triangulation<1>::substruct_iterator
+Triangulation<1>::face_iterator
CellAccessor<1>::face (const unsigned int) const {
Assert (false, ExcNotUsefulForThisDimension());
return 0;
-Triangulation<2>::substruct_iterator
+Triangulation<2>::face_iterator
CellAccessor<2>::face (const unsigned int i) const {
return line(i);
};
#include "../../../mia/vectormemory.h"
#include "../../../mia/cg.h"
+#include <map.h>
#include <algo.h>
dVector &solution,
dVector &right_hand_side,
const DirichletBC &dirichlet_bc) {
-// l;
+ 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
+ // with two bc (e.g. a corner node, with
+ // 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);
+
+ map<int,double>::const_iterator dof, endd;
+ const unsigned int n_dofs = (unsigned int)matrix.m();
+ const dSMatrixStruct &sparsity = matrix.get_sparsity_pattern();
+ const int *sparsity_rowstart = sparsity.get_rowstart_indices(),
+ *sparsity_colnums = sparsity.get_column_numbers();
+
+ for (dof=boundary_values.begin(), endd=boundary_values.end(); dof != endd; ++dof)
+ {
+ // for each boundary dof:
+
+ // set entries of this line
+ // to zero
+ for (int j=sparsity_rowstart[(*dof).first];
+ j<sparsity_rowstart[(*dof).first+1]; ++j)
+ if (sparsity_colnums[j] != (*dof).first)
+ // if not main diagonal entry
+ matrix.global_entry(j) = 0.;
+
+ // set right hand side to
+ // wanted value: if main diagonal
+ // entry nonzero, don't touch it
+ // and scale rhs accordingly. If
+ // zero, take the first main
+ // diagonal entry we can find, or
+ // one if no nonzero main diagonal
+ // element exists.
+ if (matrix.diag_element((*dof).first) != 0.0)
+ right_hand_side((*dof).first) = (*dof).second *
+ matrix.diag_element((*dof).first);
+ else
+ {
+ double first_diagonal_entry = 1;
+ for (unsigned int i=0; i<n_dofs; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ first_diagonal_entry = matrix.diag_element(i);
+ break;
+ };
+
+ matrix.set((*dof).first, (*dof).first,
+ first_diagonal_entry);
+ 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 (int j=sparsity_rowstart[row];
+ j<sparsity_rowstart[row+1]; ++row)
+ if ((sparsity_colnums[j] == (signed int)(*dof).first) &&
+ ((signed int)row != (*dof).first))
+ // this line has an entry
+ // in the regarding column
+ // but this is not the main
+ // diagonal entry
+ {
+ // correct right hand side
+ right_hand_side(row) -= matrix.global_entry(j)/diagonal_entry *
+ (*dof).second;
+
+ // set matrix entry to zero
+ matrix.global_entry(j) = 0.;
+ };
+
+
+ // preset solution vector
+ solution((*dof).first) = (*dof).second;
+ };
+};
+
+
+
+
+
+void
+ProblemBase<1>::make_boundary_value_list (const DirichletBC &,
+ map<int,double> &) const {
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+
+template <int dim>
+void
+ProblemBase<dim>::make_boundary_value_list (const DirichletBC &dirichlet_bc,
+ map<int,double> &boundary_values) const {
+ DoFHandler<dim>::active_face_iterator face = dof_handler->begin_active_face(),
+ endf = dof_handler->end_face();
+ DirichletBC::const_iterator function_ptr;
+ for (; face!=endf; ++face)
+ if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
+ dirichlet_bc.end())
+ // face is subject to one of the
+ // bc listed in #dirichlet_bc#
+ {
+ // 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->get_dof_indices (face_dofs);
+//physical location
+ (*function_ptr).second->value_list (dof_locations, dof_values);
+
+ // enter into list
+ for (unsigned int i=0; i<face_dofs.size(); ++i)
+ boundary_values[face_dofs[i]] = dof_values[i];
+ };
};
#include <grid/tria_boundary.h>
#include <grid/dof_constraints.h>
#include <basic/data_io.h>
+#include <basic/function.h>
#include <fe/fe_lib.h>
#include <fe/quadrature_lib.h>
#include <numerics/base.h>
#include <numerics/assembler.h>
#include <lac/dsmatrix.h>
+
+#include <map.h>
#include <fstream.h>
#include <cmath>
extern "C" {