*
* @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:
/**
* 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
* 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
};
+
/**
* This class is an actual implementation of the functionality proposed by
* the @ref{DataOut_DoFData} class. It offers a function @p{build_patches} that
--- /dev/null
+//---------------------------- 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
#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),
-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());
-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 ());
};
-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 ();
-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;
-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());
};
// 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
-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;
};
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
--- /dev/null
+//---------------------------- 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
+
projects / dealii-svn.git / commitdiff