QTrapez ();
};
-
-
+/**
+ * Iterated trapezoidal rule. The original reason for implementing
+ * this was gnuplot output for patches. Before this comment is
+ * removed, please do not rely on the correctness of the quadrature
+ * weights.
+ */
+template <int dim>
+class QIteratedTrapez :
+ public Quadrature<dim>
+{
+public:
+ QIteratedTrapez(const unsigned intervals);
+};
/*---------------------------- quadrature_lib.h ---------------------------*/
/* end of #ifndef __quadrature_lib_H */
};
};
+template <>
+QIteratedTrapez<1>::QIteratedTrapez (const unsigned n) :
+ Quadrature<1> (n+1)
+{
+ // Loop over INTERIOR points
+ for (unsigned int i=1; i<n; ++i)
+ {
+ quadrature_points[i] = Point<1>(1.*i/n);
+ weights[i] = 1./n;
+ };
+ quadrature_points[0] = Point<1>(0.);
+ weights[0] = .5/n;
+ quadrature_points[n] = Point<1>(1.);
+ weights[n] = .5/n;
+}
+
#endif
Assert (dim<3, ExcInternalError());
};
+template <int dim>
+QIteratedTrapez<dim>::QIteratedTrapez (const unsigned n) :
+ Quadrature<dim> (QIteratedTrapez<dim-1>(n), QIteratedTrapez<1>(n))
+{
+ Assert (dim<3, ExcInternalError());
+};
template class QMidpoint<deal_II_dimension>;
template class QSimpson<deal_II_dimension>;
template class QTrapez<deal_II_dimension>;
+template class QIteratedTrapez<deal_II_dimension>;
#endif
template <int dim, class Accessor> class TriaActiveIterator;
template <int dim> class Triangulation;
-template <int dim> class FiniteElementBase;
+template <int dim> class FiniteElement;
template <int dim> class Function;
class dVector;
* don't need them after calling this
* function, or if so store them.
*/
- virtual void distribute_dofs (const FiniteElementBase<dim> &);
+ virtual void distribute_dofs (const FiniteElement<dim> &);
/**
* Renumber the degrees of freedom according
* function is inline, so it should
* be reasonably fast.
*/
- const FiniteElementBase<dim> & get_fe () const;
+ const FiniteElement<dim> & get_fe () const;
/**
* Return a constant reference to the
* dofs per line, but also restriction
* and prolongation matrices, etc.
*/
- SmartPointer<const FiniteElementBase<dim> > selected_fe;
+ SmartPointer<const FiniteElement<dim> > selected_fe;
private:
/**
template <int dim>
inline
-const FiniteElementBase<dim> & DoFHandler<dim>::get_fe () const {
+const FiniteElement<dim> & DoFHandler<dim>::get_fe () const {
return *selected_fe;
};
* don't need them after calling this
* function, or if so store them.
*/
- virtual void distribute_dofs (const FiniteElementBase<dim> &);
+ virtual void distribute_dofs (const FiniteElement<dim> &);
/**
* Renumber the degrees of freedom according
* don't need them after calling this
* function, or if so store them.
*/
- virtual void distribute_dofs (const FiniteElementBase<dim> &);
+ virtual void distribute_dofs (const FiniteElement<dim> &);
/**
* Renumber the degrees of freedom according
/**
* Write data and grid in GNUPLOT format.
+ * Write a patch for each grid cell using
+ * #accuracy# subintervals per dimension.
*/
- void write_gnuplot (ostream &out) const;
+ void write_gnuplot (ostream &out, const unsigned accuracy=1) const;
+
+ /**
+ * Write data and grid in GNUPLOT format.
+ * Only the edges between cells are written.
+ */
+ void write_gnuplot_draft (ostream &out) const;
/**
* Write data of first vector and grid
template <int dim>
-void DoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &ff) {
+void DoFHandler<dim>::distribute_dofs (const FiniteElement<dim> &ff) {
Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
selected_fe = &ff;
template <int dim>
-void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
+void MGDoFHandler<dim>::distribute_dofs (const FiniteElement<dim> &fe) {
// first distribute global dofs
DoFHandler<dim>::distribute_dofs (fe);
template <int dim>
-void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
+void MGDoFHandler<dim>::distribute_dofs (const FiniteElement<dim> &fe) {
// first distribute global dofs
DoFHandler<dim>::distribute_dofs (fe);
#include <grid/dof.h>
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
+#include <grid/tria_boundary.h>
#include <grid/tria.h>
#include <grid/geometry_info.h>
+#include <fe/quadrature_lib.h>
+#include <fe/fe_values.h>
#include <fe/fe.h>
#include <lac/dvector.h>
template <int dim>
-void DataOut<dim>::write_gnuplot (ostream &out) const {
+void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
+{
+ Assert (dofs != 0, ExcNoDoFHandlerSelected());
+ Assert ((1<=dim) && (dim<=3), ExcNotImplemented());
+
+ if (dim==1)
+ Assert (dofs->get_fe().dofs_per_vertex==1,
+ ExcIncorrectDofsPerVertex());
+
+ // write preamble
+ if (true)
+ {
+ // block this to have local
+ // variables destroyed after
+ // use
+ time_t time1= time (0);
+ tm *time = localtime(&time1);
+ out << "# This file was generated by the deal.II library." << endl
+ << "# Date = "
+ << time->tm_year+1900 << "/"
+ << time->tm_mon+1 << "/"
+ << time->tm_mday << endl
+ << "# Time = "
+ << time->tm_hour << ":"
+ << setw(2) << time->tm_min << ":"
+ << setw(2) << time->tm_sec << endl
+ << "#" << endl
+ << "# For a description of the GNUPLOT format see the GNUPLOT manual."
+ << endl
+ << "#" << endl
+ << "# ";
+ switch (dim)
+ {
+ case 1:
+ out << "<x> ";
+ break;
+ case 2:
+ out << "<x> <y> ";
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ for (unsigned int i=0; i!=data.size(); ++i)
+ out << '<'
+ << data[i].name << ':' << data[i].units
+ << "> ";
+ out << endl;
+
+ };
+
+
+ DoFHandler<dim>::active_cell_iterator cell;
+ DoFHandler<dim>::active_cell_iterator endc = dofs->end();
+
+ QIteratedTrapez<dim> points(accuracy);
+
+ FEValues<dim> fe(dofs->get_fe(), points, UpdateFlags(update_q_points));
+ const StraightBoundary<dim> boundary;
+ vector<double> *values = new vector<double> [data.size()];
+
+ for (cell=dofs->begin_active(); cell!=endc; ++cell)
+ {
+ fe.reinit(cell, boundary);
+
+ for (unsigned i=0; i<data.size(); ++i)
+ {
+ values[i].resize(points.n_quadrature_points);
+ fe.get_function_values(*data[i].data, values[i]);
+ }
+
+ unsigned supp_pt;
+
+ switch (dim)
+ {
+ case 1:
+ // one dimension: write
+ // left vertex, right vertex
+ // and data values
+ for (supp_pt = 0; supp_pt<points.n_quadrature_points; ++supp_pt)
+ {
+ Point<dim> pt = fe.quadrature_point(supp_pt);
+ out << pt << "\t";
+ for (unsigned int i=0; i!=data.size(); ++i)
+ out << values[i][supp_pt]
+ << ' ';
+ out << endl;
+ };
+
+ break;
+
+ case 2:
+ // two dimension: output
+ // grid and data as a sequence
+ // of lines in 3d
+
+ for (unsigned xpt = 0, supp_pt = 0;
+ xpt <= accuracy; ++xpt)
+ {
+ for(unsigned ypt = 0; ypt <= accuracy; ++ypt, ++supp_pt)
+ {
+ Point<dim> pt = fe.quadrature_point(supp_pt);
+ out << pt << "\t";
+
+ for (unsigned int i=0; i!=data.size(); ++i)
+ out << values[i][supp_pt]
+ << ' ';
+ out << endl;
+ }
+ out << endl;
+ }
+
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ out << endl;
+ }
+ delete [] values;
+
+ if (!out)
+ throw GlobalExcIO ();
+}
+
+
+
+template <int dim>
+void DataOut<dim>::write_gnuplot_draft (ostream &out) const
+{
Assert (dofs != 0, ExcNoDoFHandlerSelected());
Assert ((1<=dim) && (dim<=3), ExcNotImplemented());