]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Proper FiniteElement in DofHandler
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 5 Sep 1998 15:31:08 +0000 (15:31 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 5 Sep 1998 15:31:08 +0000 (15:31 +0000)
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
deal.II/base/source/quadrature_lib.cc
deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/dofs/mg_dof_handler.h
deal.II/deal.II/include/multigrid/mg_dof_handler.h
deal.II/deal.II/include/numerics/data_io.h
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/dofs/mg_dof_handler.cc
deal.II/deal.II/source/multigrid/mg_dof_handler.cc
deal.II/deal.II/source/numerics/data_io.cc

index 1f1997150a227e493d12322e3c6bc65ba241d16a..dcdfdcb913b9188796c06ff439143bd30588a982 100644 (file)
@@ -146,8 +146,19 @@ class QTrapez : public Quadrature<dim> {
     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 */
index 597306dbb1623e0c66e96ac2f75685a6ea19678c..3e700e8a34f531c3229175e8cca874ab0440282b 100644 (file)
@@ -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<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
 
 
@@ -391,6 +407,12 @@ QTrapez<dim>::QTrapez () :
   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());
+};
 
 
 
@@ -405,6 +427,7 @@ template class QGauss8<deal_II_dimension>;
 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
index 24e2ea14e163818a4f3b72c37987bac75f66d7a9..3ed68589866d75d1eb201b425115fd8823d9b56e 100644 (file)
@@ -27,7 +27,7 @@ template <int dim, class Accessor> class TriaIterator;
 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;
@@ -498,7 +498,7 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      * 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
@@ -1242,7 +1242,7 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      * 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
@@ -1324,7 +1324,7 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      * dofs per line, but also restriction
                                      * and prolongation matrices, etc.
                                      */
-    SmartPointer<const FiniteElementBase<dim> > selected_fe;
+    SmartPointer<const FiniteElement<dim> > selected_fe;
 
   private:
                                     /**
@@ -1423,7 +1423,7 @@ class DoFHandler : public DoFDimensionInfo<dim> {
 
 template <int dim>
 inline
-const FiniteElementBase<dim> & DoFHandler<dim>::get_fe () const {
+const FiniteElement<dim> & DoFHandler<dim>::get_fe () const {
   return *selected_fe;
 };
 
index 66b22d80421d0af75804b006247f27cd628c8629..13cb03c54f23db1a4d9609725412ae3f213d9547 100644 (file)
@@ -139,7 +139,7 @@ class MGDoFHandler : public DoFHandler<dim> {
                                      * 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
index 66b22d80421d0af75804b006247f27cd628c8629..13cb03c54f23db1a4d9609725412ae3f213d9547 100644 (file)
@@ -139,7 +139,7 @@ class MGDoFHandler : public DoFHandler<dim> {
                                      * 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
index 2c28a10f9dd5484de9212604edee4e1a4d008e1d..99340e0bd7aa62c7099b7f763550bc3ca75ceb10 100644 (file)
@@ -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
index b24f9fbfad0e8744a146bced4c5a3e16ba765574..b6216e5241ee3134c53f36c2f5b50a4f378e06ac 100644 (file)
@@ -882,7 +882,7 @@ const Triangulation<dim> & DoFHandler<dim>::get_tria () const {
 
 
 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;
index 6466c1ed4277f4308d6c1838b8e5f79855e90d28..942ed376134f0615ab09a66463e4fa9e042fb598 100644 (file)
@@ -858,7 +858,7 @@ MGDoFHandler<dim>::last_active_quad () const {
 
 
 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);
 
index 6466c1ed4277f4308d6c1838b8e5f79855e90d28..942ed376134f0615ab09a66463e4fa9e042fb598 100644 (file)
@@ -858,7 +858,7 @@ MGDoFHandler<dim>::last_active_quad () const {
 
 
 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);
 
index c364746854e9a33011db57bbf22819f4647ab3e7..201916cfbcee8e335eea248c61ded2f4c306759d 100644 (file)
@@ -6,8 +6,11 @@
 #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>
 
@@ -458,7 +461,135 @@ void DataOut<dim>::write_ucd_faces (ostream &out,
 
 
 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());
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.