From be6c5dba41b192fe269917f6f66483b0380d7861 Mon Sep 17 00:00:00 2001 From: guido Date: Sat, 5 Sep 1998 15:31:08 +0000 Subject: [PATCH] Proper FiniteElement in DofHandler Great new write_gnuplot for higher order elements (old one is write_gnuplot_draft) git-svn-id: https://svn.dealii.org/trunk@576 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature_lib.h | 15 +- deal.II/base/source/quadrature_lib.cc | 23 +++ deal.II/deal.II/include/dofs/dof_handler.h | 10 +- deal.II/deal.II/include/dofs/mg_dof_handler.h | 2 +- .../include/multigrid/mg_dof_handler.h | 2 +- deal.II/deal.II/include/numerics/data_io.h | 10 +- deal.II/deal.II/source/dofs/dof_handler.cc | 2 +- deal.II/deal.II/source/dofs/mg_dof_handler.cc | 2 +- .../source/multigrid/mg_dof_handler.cc | 2 +- deal.II/deal.II/source/numerics/data_io.cc | 133 +++++++++++++++++- 10 files changed, 187 insertions(+), 14 deletions(-) diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index 1f1997150a..dcdfdcb913 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -146,8 +146,19 @@ class QTrapez : public Quadrature { 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 +class QIteratedTrapez : + public Quadrature +{ +public: + QIteratedTrapez(const unsigned intervals); +}; /*---------------------------- quadrature_lib.h ---------------------------*/ /* end of #ifndef __quadrature_lib_H */ diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 597306dbb1..3e700e8a34 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -317,6 +317,22 @@ QTrapez<1>::QTrapez () : }; }; +template <> +QIteratedTrapez<1>::QIteratedTrapez (const unsigned n) : + Quadrature<1> (n+1) +{ + // Loop over INTERIOR points + for (unsigned int i=1; i(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 @@ -391,6 +407,12 @@ QTrapez::QTrapez () : Assert (dim<3, ExcInternalError()); }; +template +QIteratedTrapez::QIteratedTrapez (const unsigned n) : + Quadrature (QIteratedTrapez(n), QIteratedTrapez<1>(n)) +{ + Assert (dim<3, ExcInternalError()); +}; @@ -405,6 +427,7 @@ template class QGauss8; template class QMidpoint; template class QSimpson; template class QTrapez; +template class QIteratedTrapez; #endif diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 24e2ea14e1..3ed6858986 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -27,7 +27,7 @@ template class TriaIterator; template class TriaActiveIterator; template class Triangulation; -template class FiniteElementBase; +template class FiniteElement; template class Function; class dVector; @@ -498,7 +498,7 @@ class DoFHandler : public DoFDimensionInfo { * don't need them after calling this * function, or if so store them. */ - virtual void distribute_dofs (const FiniteElementBase &); + virtual void distribute_dofs (const FiniteElement &); /** * Renumber the degrees of freedom according @@ -1242,7 +1242,7 @@ class DoFHandler : public DoFDimensionInfo { * function is inline, so it should * be reasonably fast. */ - const FiniteElementBase & get_fe () const; + const FiniteElement & get_fe () const; /** * Return a constant reference to the @@ -1324,7 +1324,7 @@ class DoFHandler : public DoFDimensionInfo { * dofs per line, but also restriction * and prolongation matrices, etc. */ - SmartPointer > selected_fe; + SmartPointer > selected_fe; private: /** @@ -1423,7 +1423,7 @@ class DoFHandler : public DoFDimensionInfo { template inline -const FiniteElementBase & DoFHandler::get_fe () const { +const FiniteElement & DoFHandler::get_fe () const { return *selected_fe; }; diff --git a/deal.II/deal.II/include/dofs/mg_dof_handler.h b/deal.II/deal.II/include/dofs/mg_dof_handler.h index 66b22d8042..13cb03c54f 100644 --- a/deal.II/deal.II/include/dofs/mg_dof_handler.h +++ b/deal.II/deal.II/include/dofs/mg_dof_handler.h @@ -139,7 +139,7 @@ class MGDoFHandler : public DoFHandler { * don't need them after calling this * function, or if so store them. */ - virtual void distribute_dofs (const FiniteElementBase &); + virtual void distribute_dofs (const FiniteElement &); /** * Renumber the degrees of freedom according diff --git a/deal.II/deal.II/include/multigrid/mg_dof_handler.h b/deal.II/deal.II/include/multigrid/mg_dof_handler.h index 66b22d8042..13cb03c54f 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_handler.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_handler.h @@ -139,7 +139,7 @@ class MGDoFHandler : public DoFHandler { * don't need them after calling this * function, or if so store them. */ - virtual void distribute_dofs (const FiniteElementBase &); + virtual void distribute_dofs (const FiniteElement &); /** * Renumber the degrees of freedom according diff --git a/deal.II/deal.II/include/numerics/data_io.h b/deal.II/deal.II/include/numerics/data_io.h index 2c28a10f9d..99340e0bd7 100644 --- a/deal.II/deal.II/include/numerics/data_io.h +++ b/deal.II/deal.II/include/numerics/data_io.h @@ -332,8 +332,16 @@ class DataOut { /** * 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 diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index b24f9fbfad..b6216e5241 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -882,7 +882,7 @@ const Triangulation & DoFHandler::get_tria () const { template -void DoFHandler::distribute_dofs (const FiniteElementBase &ff) { +void DoFHandler::distribute_dofs (const FiniteElement &ff) { Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); selected_fe = &ff; diff --git a/deal.II/deal.II/source/dofs/mg_dof_handler.cc b/deal.II/deal.II/source/dofs/mg_dof_handler.cc index 6466c1ed42..942ed37613 100644 --- a/deal.II/deal.II/source/dofs/mg_dof_handler.cc +++ b/deal.II/deal.II/source/dofs/mg_dof_handler.cc @@ -858,7 +858,7 @@ MGDoFHandler::last_active_quad () const { template -void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { +void MGDoFHandler::distribute_dofs (const FiniteElement &fe) { // first distribute global dofs DoFHandler::distribute_dofs (fe); diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc index 6466c1ed42..942ed37613 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc @@ -858,7 +858,7 @@ MGDoFHandler::last_active_quad () const { template -void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { +void MGDoFHandler::distribute_dofs (const FiniteElement &fe) { // first distribute global dofs DoFHandler::distribute_dofs (fe); diff --git a/deal.II/deal.II/source/numerics/data_io.cc b/deal.II/deal.II/source/numerics/data_io.cc index c364746854..201916cfbc 100644 --- a/deal.II/deal.II/source/numerics/data_io.cc +++ b/deal.II/deal.II/source/numerics/data_io.cc @@ -6,8 +6,11 @@ #include #include #include +#include #include #include +#include +#include #include #include @@ -458,7 +461,135 @@ void DataOut::write_ucd_faces (ostream &out, template -void DataOut::write_gnuplot (ostream &out) const { +void DataOut::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 << " "; + break; + case 2: + out << " "; + break; + default: + Assert (false, ExcNotImplemented()); + }; + for (unsigned int i=0; i!=data.size(); ++i) + out << '<' + << data[i].name << ':' << data[i].units + << "> "; + out << endl; + + }; + + + DoFHandler::active_cell_iterator cell; + DoFHandler::active_cell_iterator endc = dofs->end(); + + QIteratedTrapez points(accuracy); + + FEValues fe(dofs->get_fe(), points, UpdateFlags(update_q_points)); + const StraightBoundary boundary; + vector *values = new vector [data.size()]; + + for (cell=dofs->begin_active(); cell!=endc; ++cell) + { + fe.reinit(cell, boundary); + + for (unsigned i=0; i 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 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 +void DataOut::write_gnuplot_draft (ostream &out) const +{ Assert (dofs != 0, ExcNoDoFHandlerSelected()); Assert ((1<=dim) && (dim<=3), ExcNotImplemented()); -- 2.39.5