]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement a way to output faces instead of cells.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 7 Sep 2000 16:23:51 +0000 (16:23 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 7 Sep 2000 16:23:51 +0000 (16:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@3305 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/data_out.h
deal.II/deal.II/include/numerics/data_out_faces.h [new file with mode: 0644]
deal.II/deal.II/source/numerics/data_out.cc
deal.II/deal.II/source/numerics/data_out_faces.cc [new file with mode: 0644]

index e6349405d0eb9debc268d7b3cdfcc73d688bce12..fef249b50a3ac551e63182a7f5f28b9af5c57359 100644 (file)
@@ -133,8 +133,8 @@ template <int dim> class DoFHandler;
  *
  * @author Wolfgang Bangerth, 1999
  */
-template <int dof_handler_dim, int patch_dim>
-class DataOut_DoFData : public DataOutInterface<patch_dim>
+template <int dof_handler_dim, int patch_dim, int patch_space_dim=patch_dim>
+class DataOut_DoFData : public DataOutInterface<patch_dim,patch_space_dim>
 {
   public:
                                     /**
@@ -377,7 +377,7 @@ class DataOut_DoFData : public DataOutInterface<patch_dim>
                                      * in the output routines of the base
                                      * classes.
                                      */
-    vector<DataOutBase::Patch<patch_dim> > patches;
+    vector<DataOutBase::Patch<patch_dim,patch_space_dim> > patches;
 
                                     /**
                                      * Function by which the base
@@ -385,7 +385,8 @@ class DataOut_DoFData : public DataOutInterface<patch_dim>
                                      * what patches they shall write
                                      * to a file.
                                      */
-    virtual const vector<DataOutBase::Patch<patch_dim> > & get_patches () const;
+    virtual const vector<DataOutBase::Patch<patch_dim,patch_space_dim> > &
+    get_patches () const;
 
                                     /**
                                      * Virtual function through
@@ -397,6 +398,7 @@ class DataOut_DoFData : public DataOutInterface<patch_dim>
 };
 
 
+
 /**
  * This class is an actual implementation of the functionality proposed by
  * the @ref{DataOut_DoFData} class. It offers a function @p{build_patches} that
diff --git a/deal.II/deal.II/include/numerics/data_out_faces.h b/deal.II/deal.II/include/numerics/data_out_faces.h
new file mode 100644 (file)
index 0000000..16126c4
--- /dev/null
@@ -0,0 +1,160 @@
+//----------------------------  data_out_faces.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  data_out_faces.h  ---------------------------
+#ifndef __deal2__data_out_faces_h
+#define __deal2__data_out_faces_h
+
+
+#include <numerics/data_out.h>
+
+#include <string>
+#include <vector>
+
+template <int dim> class DoFHandler;
+
+
+/**
+ *
+ * @sect3{Interface}
+ *
+ * The interface of this class is copied from the @ref{DataOut}
+ * class. Furthermore, they share the common parent class
+ * @ref{DataOut_DoFData}. See the reference of these two classes for a
+ * discussion of the interface and how to extend it by deriving
+ * further classes from this class.
+ *
+ *
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <int dim>
+class DataOutFaces : public DataOut_DoFData<dim,dim-1,dim>
+{
+  public:
+                                    /**
+                                     * This is the central function
+                                     * of this class since it builds
+                                     * the list of patches to be
+                                     * written by the low-level
+                                     * functions of the base
+                                     * class. See the general
+                                     * documentation of this class
+                                     * for further information.
+                                     *
+                                     * The function supports
+                                     * multithreading, if deal.II is
+                                     * compiled in multithreading
+                                     * mode. The default number of
+                                     * threads to be used to build
+                                     * the patches is set to
+                                     * @p{multithread_info.n_default_threads}.
+                                     */
+    virtual void
+    build_patches (const unsigned int n_subdivisions = 1,
+                  const unsigned int n_threads      = multithread_info.n_default_threads);
+
+                                    /**
+                                     * Declare a way to describe a
+                                     * face which we would like to
+                                     * generate output for. The usual
+                                     * way would, of course, be to
+                                     * use an object of type
+                                     * @ref{DoFHandler}@p{<dim>::fec_iterator},
+                                     * but since we have to describe
+                                     * faces to objects of type
+                                     * @ref{FEValues}, we can only
+                                     * represent faces by pairs of a
+                                     * cell and the number of the
+                                     * face. This pair is here
+                                     * aliased to a name that is
+                                     * better to type.
+                                     */
+    typedef pair<typename DoFHandler<dim>::cell_iterator,unsigned int> FaceDescriptor;
+    
+    
+                                    /**
+                                     * Return the first face which we
+                                     * want output for. The default
+                                     * implementation returns the
+                                     * first active face on the
+                                     * boundary, but you might want
+                                     * to return another face in a
+                                     * derived class.
+                                     */
+    virtual FaceDescriptor first_face ();
+    
+                                    /**
+                                     * Return the next face after
+                                     * @p{face} which we want output
+                                     * for.  If there are no more
+                                     * face, @p{dofs->end()} shall be
+                                     * returned as the first
+                                     * component of the return value.
+                                     *
+                                     * The default implementation
+                                     * returns the next active face
+                                     * on the boundary, but you might
+                                     * want to return other faces in
+                                     * a derived class. Note that the
+                                     * default implementation assumes
+                                     * that the given @p{face} is
+                                     * active, which is guaranteed as
+                                     * long as @p{first_face} is also
+                                     * used from the default
+                                     * implementation. Overloading
+                                     * only one of the two functions
+                                     * might not be a good idea.
+                                     */
+    virtual FaceDescriptor next_face (const FaceDescriptor &face);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException1 (ExcInvalidNumberOfSubdivisions,
+                   int,
+                   << "The number of subdivisions per patch, " << arg1
+                   << ", is not valid.");
+    
+  private:
+                                    /**
+                                     * All data needed in one thread
+                                     * is gathered in the struct
+                                     * Data.
+                                     * The data is handled globally
+                                     * to avoid allocation of memory
+                                     * in the threads.
+                                     */
+    struct Data 
+    {
+       unsigned int n_threads;
+       unsigned int this_thread;
+       unsigned int n_components;
+       unsigned int n_datasets;
+       unsigned int n_subdivisions;
+       vector<double>          patch_values;
+       vector<Vector<double> > patch_values_system;
+       Data ()
+         {}
+    };
+                                    /**
+                                     * Builds every @p{n_threads}'s
+                                     * patch. This function may be
+                                     * called in parallel.
+                                     * If multithreading is not
+                                     * used, the function is called
+                                     * once and generates all patches.
+                                     */
+    void build_some_patches (Data data);
+};
+
+
+#endif
index 053bd3a70c68d78d714b9164f47ea77a9c3c1949..6db618835d9a1ab2ff2f4a8a52d73acc883817eb 100644 (file)
@@ -28,8 +28,8 @@
 
 #include <strstream>
 
-template <int dof_handler_dim, int patch_dim>
-DataOut_DoFData<dof_handler_dim,patch_dim>::DataEntry::
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::DataEntry::
 DataEntry (const Vector<double> *data,
           const vector<string> &names) :
                data(data),
@@ -38,24 +38,24 @@ DataEntry (const Vector<double> *data,
 
 
 
-template <int dof_handler_dim, int patch_dim>
-DataOut_DoFData<dof_handler_dim,patch_dim>::DataOut_DoFData () :
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::DataOut_DoFData () :
                dofs(0)
 {};
 
 
 
-template <int dof_handler_dim, int patch_dim>
-DataOut_DoFData<dof_handler_dim,patch_dim>::~DataOut_DoFData ()
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::~DataOut_DoFData ()
 {
   clear ();
 };
 
 
 
-template <int dof_handler_dim, int patch_dim>
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
 void
-DataOut_DoFData<dof_handler_dim,patch_dim>::
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::
 attach_dof_handler (const DoFHandler<dof_handler_dim> &d)
 {
   Assert (dof_data.size() == 0, ExcOldDataStillPresent());
@@ -71,9 +71,9 @@ attach_dof_handler (const DoFHandler<dof_handler_dim> &d)
 
 
 
-template <int dof_handler_dim, int patch_dim>
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
 void
-DataOut_DoFData<dof_handler_dim,patch_dim>::add_data_vector (const Vector<double> &vec,
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::add_data_vector (const Vector<double> &vec,
                                                             const vector<string> &names)
 {
   Assert (dofs != 0, ExcNoDoFHandlerSelected ());
@@ -112,9 +112,9 @@ DataOut_DoFData<dof_handler_dim,patch_dim>::add_data_vector (const Vector<double
 };
 
 
-template <int dof_handler_dim, int patch_dim>
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
 void
-DataOut_DoFData<dof_handler_dim,patch_dim>::add_data_vector (const Vector<double> &vec,
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::add_data_vector (const Vector<double> &vec,
                                                             const string         &name)
 {
   unsigned int n_components = dofs->get_fe().n_components ();
@@ -146,22 +146,22 @@ DataOut_DoFData<dof_handler_dim,patch_dim>::add_data_vector (const Vector<double
 
 
 
-template <int dof_handler_dim, int patch_dim>
-void DataOut_DoFData<dof_handler_dim,patch_dim>::clear_data_vectors ()
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
+void DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::clear_data_vectors ()
 {
   dof_data.erase (dof_data.begin(), dof_data.end());
   cell_data.erase (cell_data.begin(), cell_data.end());
 
                                   // delete patches
-  vector<DataOutBase::Patch<patch_dim> > dummy;
+  vector<DataOutBase::Patch<patch_dim,patch_space_dim> > dummy;
   patches.swap (dummy);
 }
 
 
 
-template <int dof_handler_dim, int patch_dim>
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
 void
-DataOut_DoFData<dof_handler_dim,patch_dim>::clear_input_data_references ()
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::clear_input_data_references ()
 {
   for (unsigned int i=0; i<dof_data.size(); ++i)
     dof_data[i].data = 0;
@@ -178,9 +178,9 @@ DataOut_DoFData<dof_handler_dim,patch_dim>::clear_input_data_references ()
 
 
 
-template <int dof_handler_dim, int patch_dim>
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
 void
-DataOut_DoFData<dof_handler_dim,patch_dim>::clear ()
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::clear ()
 {
   dof_data.erase (dof_data.begin(), dof_data.end());
   cell_data.erase (cell_data.begin(), cell_data.end());
@@ -192,14 +192,14 @@ DataOut_DoFData<dof_handler_dim,patch_dim>::clear ()
     };
 
                                   // delete patches
-  vector<DataOutBase::Patch<patch_dim> > dummy;
+  vector<DataOutBase::Patch<patch_dim,patch_space_dim> > dummy;
   patches.swap (dummy);
 }
 
 
-template <int dof_handler_dim, int patch_dim>
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
 vector<string>
-DataOut_DoFData<dof_handler_dim,patch_dim>::get_dataset_names () const 
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::get_dataset_names () const 
 {
   vector<string> names;
                                   // collect the names of dof
@@ -220,9 +220,9 @@ DataOut_DoFData<dof_handler_dim,patch_dim>::get_dataset_names () const
 
 
 
-template <int dof_handler_dim, int patch_dim>
-const vector<typename DataOutBase::Patch<patch_dim> > &
-DataOut_DoFData<dof_handler_dim,patch_dim>::get_patches () const
+template <int dof_handler_dim, int patch_dim, int patch_space_dim>
+const vector<typename DataOutBase::Patch<patch_dim, patch_space_dim> > &
+DataOut_DoFData<dof_handler_dim,patch_dim,patch_space_dim>::get_patches () const
 {
   return patches;
 };
@@ -430,3 +430,7 @@ template class DataOut_DoFData<deal_II_dimension,deal_II_dimension+1>;
 template class DataOut<deal_II_dimension>;
 
 
+// for 3d, also generate an extra class
+#if deal_II_dimension >= 3
+template class DataOut_DoFData<deal_II_dimension,deal_II_dimension-1,deal_II_dimension>;
+#endif
diff --git a/deal.II/deal.II/source/numerics/data_out_faces.cc b/deal.II/deal.II/source/numerics/data_out_faces.cc
new file mode 100644 (file)
index 0000000..c04b828
--- /dev/null
@@ -0,0 +1,289 @@
+//----------------------------  data_out_faces.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  data_out_faces.cc  ---------------------------
+
+
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <numerics/data_out_faces.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <fe/fe.h>
+#include <fe/fe_values.h>
+
+#ifdef DEAL_II_USE_MT
+#include <base/thread_management.h>
+#endif
+
+#include <strstream>
+
+
+
+template <int dim>
+void DataOutFaces<dim>::build_some_patches (Data data)
+{
+  QTrapez<1>        q_trapez;
+  QIterated<dim-1>  patch_points (q_trapez, data.n_subdivisions);
+  
+  FEFaceValues<dim> fe_patch_values(dofs->get_fe(),
+                                   patch_points,
+                                   update_values);
+
+  const unsigned int n_q_points = patch_points.n_quadrature_points;
+  
+  unsigned int face_number = 0;
+  typename vector<DataOutBase::Patch<dim-1,dim> >::iterator patch = patches.begin();
+  FaceDescriptor face=first_face();
+
+                                  // get first face in this thread
+  for (unsigned int i=0; (i<data.this_thread)&&(face.first != dofs->end()); ++i)
+    {
+      ++patch;
+      ++face_number;
+      face=next_face(face);
+    }
+
+                                  // now loop over all cells and
+                                  // actually create the patches
+  for (; face.first != dofs->end();)
+    {
+      Assert (patch != patches.end(), ExcInternalError());
+      
+      for (unsigned int vertex=0; vertex<GeometryInfo<dim-1>::vertices_per_cell; ++vertex)
+       patch->vertices[vertex] = face.first->face(face.second)->vertex(vertex);
+      
+      if (data.n_datasets > 0)
+       {
+         fe_patch_values.reinit (face.first, face.second);
+         
+                                          // first fill dof_data
+         for (unsigned int dataset=0; dataset<dof_data.size(); ++dataset)
+           {
+             if (data.n_components == 1)
+               {
+                 fe_patch_values.get_function_values (*dof_data[dataset].data,
+                                                      data.patch_values);
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   patch->data(dataset,q) = data.patch_values[q];
+               }
+             else
+                                                // system of components
+               {
+                 fe_patch_values.get_function_values (*dof_data[dataset].data,
+                                                      data.patch_values_system);
+                 for (unsigned int component=0; component<data.n_components;
+                      ++component)
+                   for (unsigned int q=0; q<n_q_points; ++q)
+                     patch->data(dataset*data.n_components+component,q) =
+                       data.patch_values_system[q](component);
+               };
+           };
+
+                                          // then do the cell data
+         for (unsigned int dataset=0; dataset<cell_data.size(); ++dataset)
+           {
+             const double value = (*cell_data[dataset].data)(face_number);
+             for (unsigned int q=0; q<n_q_points; ++q)
+               patch->data(dataset+dof_data.size()*data.n_components,q) =
+                 value;
+           };
+       };
+                                      // next cell (patch) in this thread
+      for (unsigned int i=0;
+          (i<data.n_threads)&&(face.first != dofs->end()); ++i)
+       {
+         ++patch;
+         ++face_number;
+         face=next_face(face);
+       }
+    };
+};
+
+
+
+
+template <int dim>
+void DataOutFaces<dim>::build_patches (const unsigned int n_subdivisions,
+                                      const unsigned int n_threads_) 
+{
+  Assert (n_subdivisions >= 1,
+         ExcInvalidNumberOfSubdivisions(n_subdivisions));
+
+  typedef DataOut_DoFData<dim,dim+1> BaseClass;
+  Assert (dofs != 0, typename BaseClass::ExcNoDoFHandlerSelected());
+
+#ifdef DEAL_II_USE_MT
+  const unsigned int n_threads = n_threads_;
+#else
+                                  // access this variable to avoid
+                                  // compiler warning about unused
+                                  // var:
+  (void)n_threads_;
+  const unsigned int n_threads = 1;
+#endif
+
+
+                                  // before we start the loop:
+                                  // create a quadrature rule that
+                                  // actually has the points on this
+                                  // patch
+  QTrapez<1>       q_trapez;
+  QIterated<dim-1> patch_points (q_trapez, n_subdivisions);
+
+  const unsigned int n_q_points     = patch_points.n_quadrature_points;
+  const unsigned int n_components   = dofs->get_fe().n_components();
+  const unsigned int n_datasets     = dof_data.size() * n_components +
+                                     cell_data.size();
+  
+                                  // clear the patches array
+  if (true)
+    {
+      vector<DataOutBase::Patch<dim-1,dim> > dummy;
+      patches.swap (dummy);
+    };
+  
+                                  // first count the cells we want to
+                                  // create patches of and make sure
+                                  // there is enough memory for that
+  unsigned int n_patches = 0;
+  for (FaceDescriptor face=first_face();
+       face.first != dofs->end();
+       face = next_face(face))
+    ++n_patches;
+
+  vector<Data> thread_data(n_threads);
+
+                                  // init data for the threads
+  for (unsigned int i=0;i<n_threads;++i)
+    {
+      thread_data[i].n_threads            = n_threads;
+      thread_data[i].this_thread          = i;
+      thread_data[i].n_components         = n_components;
+      thread_data[i].n_datasets           = n_datasets;
+      thread_data[i].n_subdivisions       = n_subdivisions;
+      thread_data[i].patch_values.resize (n_q_points);
+      thread_data[i].patch_values_system.resize (n_q_points);
+      
+      for (unsigned int k=0; k<n_q_points; ++k)
+       thread_data[i].patch_values_system[k].reinit(n_components);
+    }
+
+                                  // create the patches with default
+                                  // values. note that the evaluation
+                                  // points on the face have to be
+                                  // repeated in angular direction
+  DataOutBase::Patch<dim-1,dim>  default_patch;
+  default_patch.n_subdivisions = n_subdivisions;
+  default_patch.data.reinit (n_datasets, n_q_points);
+  patches.insert (patches.end(), n_patches, default_patch);
+
+#ifdef DEAL_II_USE_MT
+
+  Threads::ThreadManager thread_manager;  
+  for (unsigned int l=0;l<n_threads;++l)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate (&DataOutFaces<dim>::build_some_patches)
+                   .collect_args (this, thread_data[l]));
+  thread_manager.wait();
+  
+                                  // just one thread
+#else
+  build_some_patches(thread_data[0]);
+#endif
+};
+
+
+
+template <int dim>
+typename DataOutFaces<dim>::FaceDescriptor
+DataOutFaces<dim>::first_face () 
+{
+                                  // simply find first active cell
+                                  // with a face on the boundary
+  typename DoFHandler<dim>::active_cell_iterator cell = dofs->begin_active();
+  for (; cell != dofs->end(); ++cell)
+    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+      if (cell->face(f)->at_boundary())
+       return make_pair (static_cast<typename FaceDescriptor::first_type>(cell), f);
+
+                                  // ups, triangulation has no
+                                  // boundary? impossible!
+  Assert (false, ExcInternalError());
+  
+  return FaceDescriptor();
+};
+
+
+
+template <int dim>
+typename DataOutFaces<dim>::FaceDescriptor
+DataOutFaces<dim>::next_face (const FaceDescriptor &old_face)
+{
+  FaceDescriptor face = old_face;
+  
+                                  // first check whether the present
+                                  // cell has more faces on the
+                                  // boundary
+  for (unsigned int f=face.second+1; f<GeometryInfo<dim>::faces_per_cell; ++f)
+    if (face.first->face(f)->at_boundary())
+                                      // yup, that is so, so return it
+      {
+       face.second = f;
+       return face;
+      };
+
+                                  // otherwise find the next active
+                                  // cell that has a face on the
+                                  // boundary
+  
+                                  // convert the iterator to an
+                                  // active_iterator and advance
+                                  // this to the next active cell
+  typename DoFHandler<dim>::active_cell_iterator active_cell = face.first;
+
+                                  // increase face pointer by one
+  ++active_cell;
+
+                                  // while there are active cells
+  while (active_cell != dofs->end())
+    {
+                                      // check all the faces of this
+                                      // active cell
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+       if (active_cell->face(f)->at_boundary())
+         {
+           face.first  = active_cell;
+           face.second = f;
+           return face;
+         };
+                                      // the present cell had no
+                                      // faces on the boundary, so
+                                      // check next cell
+      ++active_cell;
+    };
+
+                                  // we fell off the edge, so return
+                                  // with invalid pointer
+  face.first  = dofs->end();
+  face.second = 0;
+  return face;
+};
+
+
+// explicit instantiations
+// don't instantiate anything for the 1d and 2d cases
+#if deal_II_dimension >=3
+template class DataOutFaces<deal_II_dimension>;
+#endif
+

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.