/*---------------------------- function.h ---------------------------*/
-/* <Id:> */
+/* $Id$ */
#ifndef __function_H
#define __function_H
/*---------------------------- function.h ---------------------------*/
// public classes; to be declared below
-class ParameterHandler;
-template <class UserClass> class MultipleParameterLoop;
+class ParameterHandler;
+class MultipleParameterLoop;
#include <map.h>
*/
ParameterHandler ();
+ /**
+ * Destructor. Declare this only to have
+ * a virtual destructor, which is safer
+ * as we have virtual functions.
+ * It actually does nothing spectacular.
+ */
+ virtual ~ParameterHandler ();
+
/**
* Read input from a stream until stream
* returns #eof# condition or error.
*/
const Section* get_present_changed_subsection () const;
- friend class MultipleParameterLoop;
+ friend class MultipleParameterLoop;
};
*/
MultipleParameterLoop ();
+ /**
+ * Destructor. Declare this only to have
+ * a virtual destructor, which is safer
+ * as we have virtual functions.
+ * It actually does nothing spectacular.
+ */
+ virtual ~MultipleParameterLoop ();
+
/**
* Read input from a stream until stream
* returns #eof# condition or error.
status(true) {};
+ParameterHandler::~ParameterHandler () {};
+
+
+
bool ParameterHandler::read_input (istream &input) {
String line;
int lineno=0;
ParameterHandler::Section* ParameterHandler::get_present_defaults_subsection () {
Section* sec = &defaults;
- vector<String>::const_iterator SecName;
- SecName = subsection_path.begin();
+ vector<String>::const_iterator SecName = subsection_path.begin();
while (SecName != subsection_path.end())
{
const ParameterHandler::Section* ParameterHandler::get_present_defaults_subsection () const {
- // simply call the non-const version
- // of this function
- typedef Section* (ParameterHandler::*xptr) ();
- xptr x = &get_present_defaults_subsection;
+ Section* sec = (Section*)&defaults; // not nice, but needs to be and
+ // after all: we do not change #sec#
+ vector<String>::const_iterator SecName = subsection_path.begin();
+
+ while (SecName != subsection_path.end())
+ {
+ sec = sec->subsections[*SecName];
+ ++SecName;
+ };
- return x();
+ return sec;
};
ParameterHandler::Section* ParameterHandler::get_present_changed_subsection () {
Section* sec = &changed_entries;
- vector<String>::iterator SecName;
- SecName = subsection_path.begin();
+ vector<String>::iterator SecName = subsection_path.begin();
while (SecName != subsection_path.end())
{
const ParameterHandler::Section* ParameterHandler::get_present_changed_subsection () const {
- // simply call the non-const version
- // of this function
- typedef Section* (ParameterHandler::*xptr) ();
- xptr x = &get_present_changed_subsection;
+ Section* sec = (Section*)&changed_entries; // same as in get_present_default_s...
+ vector<String>::const_iterator SecName = subsection_path.begin();
+
+ while (SecName != subsection_path.end())
+ {
+ sec = sec->subsections[*SecName];
+ ++SecName;
+ };
- return x();
+ return sec;
};
n_branches(0) {};
+MultipleParameterLoop::~MultipleParameterLoop () {};
+
+
bool MultipleParameterLoop::read_input (istream &input) {
bool x = ParameterHandler::read_input (input);
if (x) init_branches ();
LIBS = $(LIBS.g:.g=)
ifeq ($(shell uname),SunOS)
-CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
-CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
+#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \
+ -I/usr/local/source/libstdc++-2.8.0/libstdc++/stl \
+ -I/usr/local/source/libstdc++-2.8.0/libio
+CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8)
endif
Let ParameterHandler and DataIn/Out throw exceptions. Make more tests
on the input.
+Replace function objects by mem_fun, fun_ptr, ... when member function
+ templates are supported.
template <int dim> class Triangulation;
template <int dim> class DoFHandler;
-template <class T> class vector;
-
#include <grid/tria_accessor.h>
-
+#include <vector.h>
* Exception for child classes
*/
DeclException0 (ExcInvalidObject);
+ /**
+ * Exception for child classes
+ */
+ DeclException3 (ExcInvalidIndex,
+ int,
+ int,
+ int,
+ << "Invalid index " << arg1
+ << ", index must be between " << arg2
+ << " and " << arg3 << ".");
protected:
/**
template <int dim>
class DoFHandler : public DoFDimensionInfo<dim> {
public:
+ // insert these definitions for gcc2.8,
+ // since it can't inherit typedefs (I
+ // believe it should, but it can't)
+ typedef typename DoFDimensionInfo<dim>::raw_line_iterator raw_line_iterator;
+ typedef typename DoFDimensionInfo<dim>::line_iterator line_iterator;
+ typedef typename DoFDimensionInfo<dim>::active_line_iterator active_line_iterator;
+
+ typedef typename DoFDimensionInfo<dim>::raw_quad_iterator raw_quad_iterator;
+ typedef typename DoFDimensionInfo<dim>::quad_iterator quad_iterator;
+ typedef typename DoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
+
+ typedef typename DoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
+ typedef typename DoFDimensionInfo<dim>::cell_iterator cell_iterator;
+ typedef typename DoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
+
+
/**
* Constructor. Take #tria# as the
* triangulation to work on.
*/
FiniteElementBase (const FiniteElementBase &f);
-// /**
-// * Destructor. Only declared to have a
-// * virtual destructor which the compiler
-// * wants to have.
-// */
-// virtual ~FiniteElementBase () {};
+ /**
+ * Destructor. Only declared to have a
+ * virtual destructor which the compiler
+ * wants to have.
+ */
+ virtual ~FiniteElementBase () {};
/**
template <int dim>
class Triangulation : public TriaDimensionInfo<dim> {
public:
+ // insert these definitions for gcc2.8,
+ // since it can't inherit typedefs (I
+ // believe it should, but it can't)
+ typedef typename TriaDimensionInfo<dim>::raw_line_iterator raw_line_iterator;
+ typedef typename TriaDimensionInfo<dim>::line_iterator line_iterator;
+ typedef typename TriaDimensionInfo<dim>::active_line_iterator active_line_iterator;
+
+ typedef typename TriaDimensionInfo<dim>::raw_quad_iterator raw_quad_iterator;
+ typedef typename TriaDimensionInfo<dim>::quad_iterator quad_iterator;
+ typedef typename TriaDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
+
+ typedef typename TriaDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
+ typedef typename TriaDimensionInfo<dim>::cell_iterator cell_iterator;
+ typedef typename TriaDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
+
/**
* Create a triangulation and create
* the first level of the hierarchy.
@author Wolfgang Bangerth, 1998
*/
template <int dim, class Accessor>
-class TriaRawIterator : public bidirectional_iterator<Accessor,int> {
+class TriaRawIterator : public bidirectional_iterator<Accessor,int>{
public:
/**
* Empty constructor. Such an object
template <int dim> class FiniteElement;
template <int dim> class Quadrature;
template <int dim> class DataOut;
-
+template <int dim> class Function;
\item Condense the system matrix and right hand side with the constraints
induced by hanging nodes.
\end{itemize}
+
+
+ {\bf Solving}
+
+ Calling the #solve# function with a solver object, the system of equations
+ which results after having called the #assemble# function is solved. After
+ this, the solution vector is distributed again, i.e. the constrained nodes
+ are given their correct values.
*/
template <int dim>
class ProblemBase {
ProblemBase (Triangulation<dim> *tria,
DoFHandler<dim> *dof_handler);
+ /**
+ * Destructor. Declare this only to have
+ * a virtual destructor, which is safer
+ * as we have virtual functions.
+ * It actually does nothing spectacular.
+ */
+ virtual ~ProblemBase ();
+
/**
* Initiate the process of assemblage of
* vectors and system matrix. Use the
const FiniteElement<dim> &fe);
- void solve ();
- void fill_data (DataOut<dim> &) const;
+ /**
+ * Solve the system of equations.
+ */
+ virtual void solve ();
+
+ /**
+ * Integrate the difference between
+ * the solution computed before and
+ * the reference solution, which
+ * is given as a continuous function
+ * object. The integration is
+ * performed using the given quadrature
+ * rule and assumes that the given
+ * finite element objects equals
+ * that used for the computation of
+ * the solution.
+ *
+ * The result ist stored in the
+ * #difference# object, where each
+ * entry equals the L_1 norm of
+ * the difference on one cell. The
+ * order of entries is the same as
+ * a #cell_iterator# takes when
+ * started with #begin_active# and
+ * promoted with the #++# operator.
+ *
+ * You can use the
+ * #distribute_cell_to_dof_vector#
+ * function to convert cell based
+ * data to a data vector with values
+ * on the degrees of freedom, which
+ * can then be attached to a
+ * #DataOut# object.
+ *
+ * To get the global L_1 error,
+ * you have to sum up the entries
+ * in #difference#, e.g. using the STL
+ * function
+ * #accumulate (d.begin(), d.end(), 0)#,
+ * if #d# is the difference vector.
+ */
+ void integrate_L1_difference (const Function<dim> &exact_solution,
+ vector<double> &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe) const;
+
+ /**
+ * Initialize the #DataOut# object with
+ * the grid and DoF handler used for this
+ * computation, as well as with the
+ * solution. Overload this function if
+ * you have multiple data sets to be
+ * written, or alternativelt call this
+ * function and attach the additional
+ * vectors directly to the #DataOut#
+ * object.
+ *
+ * Solution name and physical units are
+ * derived by calling the virtual
+ * function #get_solution_name#.
+ */
+ virtual void fill_data (DataOut<dim> &) const;
+
+
+ /**
+ * Return solution name and
+ * physical units as a pair of
+ * #char*#. The default implementation
+ * returns #make_pair ("solution","")#,
+ * which results in "<dimensionless>"
+ * upon output.
+ * Overload this function, if you
+ * want anything else.
+ */
+ virtual pair<char*,char*> get_solution_name () const;
/**
* Exception
* Exception
*/
DeclException0 (ExcNoMemory);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidFE);
protected:
/**
class dVector;
-template <class Key, class T, class Compare> class map;
+template <class Key, class T, class Compare, class Alloc = alloc> class map;
template <class T> struct less;
#ifdef DEBUG
if (q.state() != past_the_end)
- Assert (q->used(), ExcUnusedCellAsChild());
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
#endif
return q;
};
#ifdef DEBUG
if (q.state() != past_the_end)
- Assert (q->used(), ExcUnusedCellAsChild());
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
#endif
return q;
};
#ifdef DEBUG
if (q.state() != past_the_end)
- Assert (q->used(), ExcUnusedCellAsNeighbor());
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsNeighbor());
#endif
return q;
};
#ifdef DEBUG
if (q.state() != past_the_end)
- Assert (q->used(), ExcUnusedCellAsChild());
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
#endif
return q;
};
// otherwise, the number states which
// line in the constraint matrix handles
// this index
+ //
+ // for gcc2.9: replace this by
+ // distribute(sparsity.n_rows(), -1)
vector<int> distribute;
distribute.reserve (sparsity.n_rows());
- distribute.insert (distribute.end(), sparsity.n_rows(), -1);
-
+ // replace the following line, since
+ // gcc2.8 chokes over it.
+// distribute.insert (distribute.end(), sparsity.n_rows(), -1);
+ for (unsigned int i=0; i<(unsigned int)sparsity.n_rows(); ++i)
+ distribute.push_back (-1);
+
for (int c=0; c<(signed int)lines.size(); ++c)
distribute[lines[c].line] = c;
// otherwise, the number states which
// line in the constraint matrix handles
// this index
+ //
+ // for gcc2.9: replace this by
+ // distribute(sparsity.n_rows(), -1)
vector<int> distribute;
distribute.reserve (sparsity.n_rows());
- distribute.insert (distribute.end(), sparsity.n_rows(), -1);
+ // replace the following line, since
+ // gcc2.8 chokes over it.
+// distribute.insert (distribute.end(), sparsity.n_rows(), -1);
+ for (unsigned int i=0; i<(unsigned int)sparsity.n_rows(); ++i)
+ distribute.push_back (-1);
+
for (int c=0; c<(signed int)lines.size(); ++c)
distribute[lines[c].line] = c;
-
template <int dim>
void DoFHandler<dim>::renumber_dofs (const RenumberingMethod method,
bool use_constraints,
make_constraint_matrix (constraints);
constraints.condense (sparsity);
};
-
+
int n_dofs = sparsity.n_rows();
// store the new dof numbers; -1 means
// that no new number was chosen yet
- vector<int> new_number(sparsity.n_rows(), -1);
-
+ //
+ // the commented line is what would be the
+ // correct way to do, but gcc2.8 chokes
+ // over that. The other lines are a
+ // workaround
+// vector<int> new_number(sparsity.n_rows(), -1);
+ vector<int> new_number;
+ new_number.reserve (sparsity.n_rows());
+ for (unsigned int i=0; i<(unsigned int)sparsity.n_rows(); ++i)
+ new_number.push_back (-1);
+
// store the indices of the dofs renumbered
// in the last round. Default to starting
// points
vector<int> last_round_dofs (starting_points);
+
// delete disallowed elements
for (unsigned int i=0; i<last_round_dofs.size(); ++i)
if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
last_round_dofs[i] = -1;
+
remove_if (last_round_dofs.begin(), last_round_dofs.end(),
bind2nd(equal_to<int>(), 0));
-
+
// now if no valid points remain:
// find dof with lowest coordination
// number
+
if (last_round_dofs.size() == 0)
{
int starting_point = -1;
sparsity.add (dofs_on_this_cell[i],
dofs_on_this_cell[j]);
};
-
+
delete[] dofs_on_this_cell;
};
template <int dim>
-void DoFHandler<dim>::transfer_cell (const TriaIterator<dim,DoFCellAccessor<dim> > &old_cell,
- const TriaIterator<dim,DoFCellAccessor<dim> > &new_cell,
+void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterator &old_cell,
+ const typename DoFHandler<dim>::cell_iterator &new_cell,
dSMatrixStruct &transfer_pattern) const {
if (!new_cell->active() && !old_cell->active())
// both cells have children; go deeper
template <int dim>
-void DoFHandler<dim>::transfer_cell (const TriaIterator<dim,DoFCellAccessor<dim> > &old_cell,
- const TriaIterator<dim,DoFCellAccessor<dim> > &new_cell,
+void DoFHandler<dim>::transfer_cell (const typename DoFHandler<dim>::cell_iterator &old_cell,
+ const typename DoFHandler<dim>::cell_iterator &new_cell,
dSMatrix &transfer_matrix) const {
if (!new_cell->active() && !old_cell->active())
// both cells have children; go deeper
-funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE)
ifeq ($(shell uname),SunOS)
-CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
-CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
+#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \
+ -I/usr/local/source/libstdc++-2.8.0/libstdc++/stl \
+ -I/usr/local/source/libstdc++-2.8.0/libio
+CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8)
endif
FEQuadratic<dim>::shape_value (const unsigned int i,
const Point<dim> &) const
{
- Assert (i<total_dofs, ExcInvalidIndex(i));
- Assert (false, ExcNotImplemented());
+ Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
+ Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
return 0.;
};
FEQuadratic<dim>::shape_grad (const unsigned int i,
const Point<dim> &) const
{
- Assert (i<total_dofs, ExcInvalidIndex(i));
- Assert (false, ExcNotImplemented());
+ Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
+ Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
return Point<dim> ();
};
const vector<Point<2> > &,
vector<dFMatrix> &,
vector<Point<2> > &) const {
- Assert (false, ExcNotImplemented());
+ Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
FECubic<dim>::shape_value (const unsigned int i,
const Point<dim> &) const
{
- Assert (i<total_dofs, ExcInvalidIndex(i));
- Assert (false, ExcNotImplemented());
+ Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
+ Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
return 0.;
};
FECubic<dim>::shape_grad (const unsigned int i,
const Point<dim> &) const
{
- Assert (i<total_dofs, ExcInvalidIndex(i));
- Assert (false, ExcNotImplemented());
+ Assert (i<total_dofs, typename FiniteElementBase<dim>::ExcInvalidIndex(i));
+ Assert (false, typename FiniteElementBase<dim>::ExcNotImplemented());
return Point<dim> ();
};
const vector<Point<2> > &,
vector<dFMatrix> &,
vector<Point<2> > &) const {
- Assert (false, ExcNotImplemented());
+ Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
};
-funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE)
ifeq ($(shell uname),SunOS)
-CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
-CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
+#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \
+ -I/usr/local/source/libstdc++-2.8.0/libstdc++/stl \
+ -I/usr/local/source/libstdc++-2.8.0/libio
+CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8)
endif
vertices.push_back (Point<1> (left));
vertices.push_back (Point<1> (right));
- vertices_used.insert (vertices_used.end(), 2, true);
+ vertices_used.push_back (true);
+ vertices_used.push_back (true);
levels[0]->lines.lines.push_back (Line(0,1));
levels[0]->lines.children.push_back (-1);
levels[0]->neighbors.push_back (make_pair(-1,-1));
levels[0]->neighbors.push_back (make_pair(-1,-1));
- levels[0]->refine_flags.insert (levels[0]->refine_flags.end(), 1, false);
+ levels[0]->refine_flags.push_back (false);
};
Assert (vertices.size() == 0, ExcTriangulationNotEmpty());
Assert (n_lines() == 0, ExcTriangulationNotEmpty());
Assert (n_quads() == 0, ExcTriangulationNotEmpty());
-
- // create vertices
- vertices.push_back (Point<2> (left,left));
- vertices.push_back (Point<2> (right,left));
- vertices.push_back (Point<2> (right,right));
- vertices.push_back (Point<2> (left,right));
- vertices_used.insert (vertices_used.end(), 4, true);
-
- // create lines
- levels[0]->lines.lines.push_back (Line (0,1));
- levels[0]->lines.children.push_back (-1);
- levels[0]->lines.used.push_back (true);
- levels[0]->lines.user_flags.push_back (false);
- levels[0]->lines.lines.push_back (Line (1,2));
- levels[0]->lines.children.push_back (-1);
- levels[0]->lines.used.push_back (true);
- levels[0]->lines.user_flags.push_back (false);
-
- levels[0]->lines.lines.push_back (Line (3,2));
- levels[0]->lines.children.push_back (-1);
- levels[0]->lines.used.push_back (true);
- levels[0]->lines.user_flags.push_back (false);
-
- levels[0]->lines.lines.push_back (Line (0,3));
- levels[0]->lines.children.push_back (-1);
- levels[0]->lines.used.push_back (true);
- levels[0]->lines.user_flags.push_back (false);
+ const Point<2> vertices[4] = { Point<2>(left,left),
+ Point<2>(right,left),
+ Point<2>(right,right),
+ Point<2>(left,right) };
+ const int cell_vertices[1][4] = { { 0,1,2,3 } };
+ vector<vector<int> > cells (1, vector<int>());
+ cells[0].reserve (4);
+ for (unsigned int j=0; j<4; ++j)
+ cells[0].push_back (cell_vertices[0][j]);
- // create quads
- levels[0]->quads.quads.push_back (Quad (0,1,2,3));
- levels[0]->quads.children.push_back (-1);
- levels[0]->quads.used.push_back (true);
- levels[0]->quads.user_flags.push_back (false);
-
- levels[0]->neighbors.push_back (make_pair(-1,-1));
- levels[0]->neighbors.push_back (make_pair(-1,-1));
- levels[0]->neighbors.push_back (make_pair(-1,-1));
- levels[0]->neighbors.push_back (make_pair(-1,-1));
-
- levels[0]->refine_flags.insert (levels[0]->refine_flags.end(), 1, false);
+ create_triangulation (vector<Point<2> >(&vertices[0], &vertices[8]),
+ cells);
};
const int cell_vertices[3][4] = {{0, 1, 4, 3},
{1, 2, 5, 4},
{3, 4, 7, 6}};
- vector<vector<int> > cells (3, vector<int>(4,-1));
+
+ // with gcc2.9, uncomment this line
+ // and delete following workaraound
+// vector<vector<int> > cells (3, vector<int>(4,-1));
+ vector<int> tmp;
+ tmp.reserve (4);
+ for (unsigned int i=0; i<4; ++i)
+ tmp.push_back (-1);
+ vector<vector<int> > cells (3, tmp);
+
for (unsigned int i=0; i<3; ++i)
for (unsigned int j=0; j<4; ++j)
cells[i][j] = cell_vertices[i][j];
{2, 3, 5, 4},
{1, 7, 5, 3},
{6, 4, 5, 7}};
- vector<vector<int> > cells (5, vector<int>(4,-1));
+
+ // with gcc2.9, uncomment this line
+ // and delete following workaraound
+// vector<vector<int> > cells (5, vector<int>(4,-1));
+ vector<int> tmp;
+ tmp.reserve (4);
+ for (unsigned int i=0; i<4; ++i)
+ tmp.push_back (-1);
+ vector<vector<int> > cells (5, tmp);
+
for (unsigned int i=0; i<5; ++i)
for (unsigned int j=0; j<4; ++j)
cells[i][j] = cell_vertices[i][j];
template <int dim>
void Triangulation<dim>::load_refine_flags (istream &in) {
- int magic_number;
+ unsigned int magic_number;
in >> magic_number;
Assert (magic_number==mn_tria_refine_flags_begin, ExcGridReadError());
template <int dim>
void Triangulation<dim>::load_user_flags_line (istream &in) {
- int magic_number;
+ unsigned int magic_number;
in >> magic_number;
Assert (magic_number==mn_tria_line_user_flags_begin, ExcGridReadError());
template <int dim>
void Triangulation<dim>::load_user_flags_quad (istream &in) {
- int magic_number;
+ unsigned int magic_number;
in >> magic_number;
Assert (magic_number==mn_tria_quad_user_flags_begin, ExcGridReadError());
if (cell->neighbor(1).state() != valid)
second_child->set_neighbor (1, cell->neighbor(1));
else
- if (cell->neighbor(1)->active)
+ if (cell->neighbor(1)->active())
second_child->set_neighbor (1, cell->neighbor(1));
else
// right neighbor is refined
-/*------------------------ Functions: CellAccessor<dim> ---------------------*/
-template <int dim>
-int CellAccessor<dim>::neighbor_index (const unsigned int i) const {
- Assert (i<2*dim, ExcInvalidNeighbor(i));
- return tria->levels[present_level]->neighbors[present_index*2*dim+i].second;
+
+/*------------------------ Functions: CellAccessor<1> -----------------------*/
+
+
+bool CellAccessor<1>::at_boundary () const {
+ return at_boundary(0) || at_boundary(1);
};
-template <int dim>
-int CellAccessor<dim>::neighbor_level (const unsigned int i) const {
- Assert (i<2*dim, ExcInvalidNeighbor(i));
- return tria->levels[present_level]->neighbors[present_index*2*dim+i].first;
+
+int CellAccessor<1>::neighbor_index (const unsigned int i) const {
+ Assert (i<2*1, ExcInvalidNeighbor(i));
+ return tria->levels[present_level]->neighbors[present_index*2*1+i].second;
};
-template <int dim>
-void CellAccessor<dim>::set_neighbor (const unsigned int i,
- const TriaIterator<dim,CellAccessor<dim> > &pointer) const {
- Assert (i<2*dim, ExcInvalidNeighbor(i));
- tria->levels[present_level]->neighbors[present_index*2*dim+i].first
+
+int CellAccessor<1>::neighbor_level (const unsigned int i) const {
+ Assert (i<2*1, ExcInvalidNeighbor(i));
+ return tria->levels[present_level]->neighbors[present_index*2*1+i].first;
+};
+
+
+
+
+void CellAccessor<1>::set_neighbor (const unsigned int i,
+ const TriaIterator<1,CellAccessor<1> > &pointer) const {
+ Assert (i<2*1, ExcInvalidNeighbor(i));
+ tria->levels[present_level]->neighbors[present_index*2*1+i].first
= pointer.accessor.present_level;
- tria->levels[present_level]->neighbors[present_index*2*dim+i].second
+ tria->levels[present_level]->neighbors[present_index*2*1+i].second
= pointer.accessor.present_index;
};
-template <int dim>
-bool CellAccessor<dim>::at_boundary (const unsigned int i) const {
+
+bool CellAccessor<1>::at_boundary (const unsigned int i) const {
Assert (used(), ExcCellNotUsed());
- Assert (i<2*dim, ExcInvalidIndex (i,0,2*dim-1));
+ Assert (i<2*1, ExcInvalidIndex (i,0,2*1-1));
return (neighbor(i).state() != valid);
};
-template <int dim>
-bool CellAccessor<dim>::refine_flag_set () const {
+
+bool CellAccessor<1>::refine_flag_set () const {
Assert (used() && active(), ExcRefineCellNotActive());
return tria->levels[present_level]->refine_flags[present_index];
};
-template <int dim>
-void CellAccessor<dim>::set_refine_flag () const {
+
+void CellAccessor<1>::set_refine_flag () const {
Assert (used() && active(), ExcRefineCellNotActive());
tria->levels[present_level]->refine_flags[present_index] = true;
};
-template <int dim>
-void CellAccessor<dim>::clear_refine_flag () const {
+
+void CellAccessor<1>::clear_refine_flag () const {
Assert (used() && active(), ExcRefineCellNotActive());
tria->levels[present_level]->refine_flags[present_index] = false;
};
-template <int dim>
-TriaIterator<dim,CellAccessor<dim> >
-CellAccessor<dim>::neighbor (const unsigned int i) const {
- TriaIterator<dim,CellAccessor<dim> > q (tria, neighbor_level (i), neighbor_index (i));
+
+TriaIterator<1,CellAccessor<1> >
+CellAccessor<1>::neighbor (const unsigned int i) const {
+ TriaIterator<1,CellAccessor<1> > q (tria, neighbor_level (i), neighbor_index (i));
#ifdef DEBUG
if (q.state() != past_the_end)
-template <int dim>
-TriaIterator<dim,CellAccessor<dim> >
-CellAccessor<dim>::child (const unsigned int i) const {
- TriaIterator<dim,CellAccessor<dim> > q (tria, present_level+1, child_index (i));
+
+TriaIterator<1,CellAccessor<1> >
+CellAccessor<1>::child (const unsigned int i) const {
+ TriaIterator<1,CellAccessor<1> > q (tria, present_level+1, child_index (i));
#ifdef DEBUG
if (q.state() != past_the_end)
-template <int dim>
-bool CellAccessor<dim>::active () const {
+
+bool CellAccessor<1>::active () const {
return !has_children();
};
-/*------------------------ Functions: CellAccessor<1> -----------------------*/
+/*------------------------ Functions: CellAccessor<2> -----------------------*/
-bool CellAccessor<1>::at_boundary () const {
- return at_boundary(0) || at_boundary(1);
+bool CellAccessor<2>::at_boundary () const {
+ return at_boundary(0) || at_boundary(1) || at_boundary(2) || at_boundary(3);
};
-/*------------------------ Functions: CellAccessor<2> -----------------------*/
+int CellAccessor<2>::neighbor_index (const unsigned int i) const {
+ Assert (i<2*2, ExcInvalidNeighbor(i));
+ return tria->levels[present_level]->neighbors[present_index*2*2+i].second;
+};
-bool CellAccessor<2>::at_boundary () const {
- return at_boundary(0) || at_boundary(1) || at_boundary(2) || at_boundary(3);
+
+
+int CellAccessor<2>::neighbor_level (const unsigned int i) const {
+ Assert (i<2*2, ExcInvalidNeighbor(i));
+ return tria->levels[present_level]->neighbors[present_index*2*2+i].first;
};
+void CellAccessor<2>::set_neighbor (const unsigned int i,
+ const TriaIterator<2,CellAccessor<2> > &pointer) const {
+ Assert (i<2*2, ExcInvalidNeighbor(i));
+ tria->levels[present_level]->neighbors[present_index*2*2+i].first
+ = pointer.accessor.present_level;
+ tria->levels[present_level]->neighbors[present_index*2*2+i].second
+ = pointer.accessor.present_index;
+};
+
+
+
+
+bool CellAccessor<2>::at_boundary (const unsigned int i) const {
+ Assert (used(), ExcCellNotUsed());
+ Assert (i<2*2, ExcInvalidIndex (i,0,2*2-1));
+
+ return (neighbor(i).state() != valid);
+};
+
+
+
+
+bool CellAccessor<2>::refine_flag_set () const {
+ Assert (used() && active(), ExcRefineCellNotActive());
+ return tria->levels[present_level]->refine_flags[present_index];
+};
+
+
+
+
+void CellAccessor<2>::set_refine_flag () const {
+ Assert (used() && active(), ExcRefineCellNotActive());
+ tria->levels[present_level]->refine_flags[present_index] = true;
+};
+
+
+
+
+void CellAccessor<2>::clear_refine_flag () const {
+ Assert (used() && active(), ExcRefineCellNotActive());
+ tria->levels[present_level]->refine_flags[present_index] = false;
+};
+
+
+
+
+TriaIterator<2,CellAccessor<2> >
+CellAccessor<2>::neighbor (const unsigned int i) const {
+ TriaIterator<2,CellAccessor<2> > q (tria, neighbor_level (i), neighbor_index (i));
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), ExcUnusedCellAsNeighbor());
+#endif
+ return q;
+};
+
+
+
+
+TriaIterator<2,CellAccessor<2> >
+CellAccessor<2>::child (const unsigned int i) const {
+ TriaIterator<2,CellAccessor<2> > q (tria, present_level+1, child_index (i));
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), ExcUnusedCellAsChild());
+#endif
+ return q;
+};
+
+
+
+
+bool CellAccessor<2>::active () const {
+ return !has_children();
+};
+
-funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE)
ifeq ($(shell uname),SunOS)
-CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
-CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
+#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \
+ -I/usr/local/source/libstdc++-2.8.0/libstdc++/stl \
+ -I/usr/local/source/libstdc++-2.8.0/libio
+CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8)
endif
#include <grid/dof_constraints.h>
#include <grid/tria_iterator.h>
#include <basic/data_io.h>
+#include <basic/function.h>
#include "../../../mia/control.h"
#include "../../../mia/vectormemory.h"
+
+
extern TriaActiveIterator<1,CellAccessor<1> > __dummy1233; // for gcc2.7
extern TriaActiveIterator<2,CellAccessor<2> > __dummy1234;
+template <int dim>
+ProblemBase<dim>::~ProblemBase () {};
+
+
template <int dim>
void ProblemBase<dim>::assemble (const Equation<dim> &equation,
+/*
+template <int dim>
+void ProblemBase<dim>::integrate_L1_difference (const Function<dim> &exact_solution,
+ vector<double> &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe) const {
+ Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE());
+
+ difference.erase (difference.begin(), difference.end());
+ difference.reserve (tria->n_cells());
+
+ // loop over all cells
+ DoFHandler<dim>::active_dof_iterator cell = dof_handler->begin_active(),
+ endc = dof_handler->end();
+ for (; cell != end; ++cell)
+ {
+ double diff=0;
+
+ difference.push_back (diff);
+ };
+};
+*/
+
+
+
template <int dim>
void ProblemBase<dim>::fill_data (DataOut<dim> &out) const {
out.clear_data_vectors ();
out.attach_dof_handler (*dof_handler);
- out.add_data_vector (solution, "solution");
+
+ pair<char*,char*> solution_name = get_solution_name ();
+ out.add_data_vector (solution,
+ solution_name.first, solution_name.second);
+};
+
+
+/*
+template <int dim>
+pair<char*,char*> ProblemBase<dim>::get_solution_name () const {
+ return make_pair("solution", "");
};
+*/
LIBS = $(LIBS.g:.g=)
ifeq ($(shell uname),SunOS)
-CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
-CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
+#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
+INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \
+ -I/usr/local/source/libstdc++-2.8.0/libstdc++/stl \
+ -I/usr/local/source/libstdc++-2.8.0/libio
+CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8)
endif