* usually @p{build_patches} or the like, which fills the @p{patches} array of
* this class.
*
+ * Regarding the templates of this class, it needs two values: first
+ * the space dimension in which the triangulation and the DoF handler
+ * operate, second the space dimension of the output to be
+ * generated. Although in most cases they are equal, there are also
+ * classes for which this does not hold, for example if one outputs
+ * the result of a computation exploiting rotational symmetry in the
+ * original domain (in which the space dimension of the output would
+ * be one higher than that of the DoF handler, see the
+ * @ref{DataOut_Rotation} class), or one might conceive that one could
+ * write a class that only outputs the solution on a cut through the
+ * domain, in which case the space dimension of the output is less
+ * than that of the DoF handler.
+ *
* @author Wolfgang Bangerth, 1999
*/
-template <int dim>
-class DataOut_DoFData : public DataOutInterface<dim>
+template <int dof_handler_dim, int patch_dim>
+class DataOut_DoFData : public DataOutInterface<patch_dim>
{
public:
/**
* and the mapping between nodes
* and node values.
*/
- void attach_dof_handler (const DoFHandler<dim> &);
+ void attach_dof_handler (const DoFHandler<dof_handler_dim> &);
/**
* Add a data vector together
/**
* Pointer to the dof handler object.
*/
- const DoFHandler<dim> *dofs;
+ const DoFHandler<dof_handler_dim> *dofs;
/**
* List of data elements with vectors of
* in the output routines of the base
* classes.
*/
- vector<DataOutBase::Patch<dim> > patches;
+ vector<DataOutBase::Patch<patch_dim> > patches;
/**
* Function by which the base
* what patches they shall write
* to a file.
*/
- virtual const vector<DataOutBase::Patch<dim> > & get_patches () const;
+ virtual const vector<DataOutBase::Patch<patch_dim> > & get_patches () const;
/**
* Virtual function through
* @author Wolfgang Bangerth, 1999
*/
template <int dim>
-class DataOut : public DataOut_DoFData<dim>
+class DataOut : public DataOut_DoFData<dim,dim>
{
public:
/**
--- /dev/null
+//---------------------------- data_out_rotation.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_rotation.h ---------------------------
+#ifndef __deal2__data_out_rotation_h
+#define __deal2__data_out_rotation_h
+
+
+#include <numerics/data_out.h>
+
+#include <string>
+#include <vector>
+
+template <int dim> class DoFHandler;
+
+/**
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <int dim>
+class DataOutRotation : public DataOut_DoFData<dim,dim+1>
+{
+ 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_patches_per_circle,
+ const unsigned int n_subdivisions = 1,
+ const unsigned int n_threads = multithread_info.n_default_threads);
+
+ /**
+ * Return the first cell which we
+ * want output for. The default
+ * implementation returns the first
+ * active cell, but you might want
+ * to return other cells in a derived
+ * class.
+ */
+ virtual typename DoFHandler<dim>::cell_iterator
+ first_cell ();
+
+ /**
+ * Return the next cell after @p{cell} which
+ * we want output for.
+ * If there are no more cells,
+ * @p{dofs->end()} shall be returned.
+ *
+ * The default
+ * implementation returns the next
+ * active cell, but you might want
+ * to return other cells in a derived
+ * class. Note that the default
+
+ * implementation assumes that
+ * the given @p{cell} is active, which
+ * is guaranteed as long as @p{first_cell}
+ * is also used from the default
+ * implementation. Overloading only one
+ * of the two functions might not be
+ * a good idea.
+ */
+ virtual typename DoFHandler<dim>::cell_iterator
+ next_cell (const typename DoFHandler<dim>::cell_iterator &cell);
+
+ /**
+ * 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_patches_per_circle;
+ 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 dim>
-DataOut_DoFData<dim>::DataEntry::DataEntry (const Vector<double> *data,
- const vector<string> &names) :
+template <int dof_handler_dim, int patch_dim>
+DataOut_DoFData<dof_handler_dim,patch_dim>::DataEntry::
+DataEntry (const Vector<double> *data,
+ const vector<string> &names) :
data(data),
names(names)
{};
-template <int dim>
-DataOut_DoFData<dim>::DataOut_DoFData () :
+template <int dof_handler_dim, int patch_dim>
+DataOut_DoFData<dof_handler_dim,patch_dim>::DataOut_DoFData () :
dofs(0)
{};
-template <int dim>
-DataOut_DoFData<dim>::~DataOut_DoFData ()
+template <int dof_handler_dim, int patch_dim>
+DataOut_DoFData<dof_handler_dim,patch_dim>::~DataOut_DoFData ()
{
clear ();
};
-template <int dim>
-void DataOut_DoFData<dim>::attach_dof_handler (const DoFHandler<dim> &d)
+template <int dof_handler_dim, int patch_dim>
+void
+DataOut_DoFData<dof_handler_dim,patch_dim>::
+attach_dof_handler (const DoFHandler<dof_handler_dim> &d)
{
Assert (dof_data.size() == 0, ExcOldDataStillPresent());
Assert (cell_data.size() == 0, ExcOldDataStillPresent());
-template <int dim>
-void DataOut_DoFData<dim>::add_data_vector (const Vector<double> &vec,
- const vector<string> &names)
+template <int dof_handler_dim, int patch_dim>
+void
+DataOut_DoFData<dof_handler_dim,patch_dim>::add_data_vector (const Vector<double> &vec,
+ const vector<string> &names)
{
Assert (dofs != 0, ExcNoDoFHandlerSelected ());
};
-template <int dim>
-void DataOut_DoFData<dim>::add_data_vector (const Vector<double> &vec,
- const string &name)
+template <int dof_handler_dim, int patch_dim>
+void
+DataOut_DoFData<dof_handler_dim,patch_dim>::add_data_vector (const Vector<double> &vec,
+ const string &name)
{
unsigned int n_components = dofs->get_fe().n_components ();
-template <int dim>
-void DataOut_DoFData<dim>::clear_data_vectors ()
+template <int dof_handler_dim, int patch_dim>
+void DataOut_DoFData<dof_handler_dim,patch_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<dim> > dummy;
+ vector<DataOutBase::Patch<patch_dim> > dummy;
patches.swap (dummy);
}
-template <int dim>
-void DataOut_DoFData<dim>::clear_input_data_references ()
+template <int dof_handler_dim, int patch_dim>
+void
+DataOut_DoFData<dof_handler_dim,patch_dim>::clear_input_data_references ()
{
for (unsigned int i=0; i<dof_data.size(); ++i)
dof_data[i].data = 0;
-
-template <int dim>
-void DataOut_DoFData<dim>::clear ()
+template <int dof_handler_dim, int patch_dim>
+void
+DataOut_DoFData<dof_handler_dim,patch_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<dim> > dummy;
+ vector<DataOutBase::Patch<patch_dim> > dummy;
patches.swap (dummy);
}
-template <int dim>
-vector<string> DataOut_DoFData<dim>::get_dataset_names () const
+template <int dof_handler_dim, int patch_dim>
+vector<string>
+DataOut_DoFData<dof_handler_dim,patch_dim>::get_dataset_names () const
{
vector<string> names;
// collect the names of dof
-template <int dim>
-const vector<typename DataOutBase::Patch<dim> > &
-DataOut_DoFData<dim>::get_patches () 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
{
return patches;
};
+/* ---------------------------------------------------------------------- */
+
+
+
template <int dim>
void DataOut<dim>::build_some_patches (Data data)
{
{
Assert (n_subdivisions >= 1,
ExcInvalidNumberOfSubdivisions(n_subdivisions));
-
- Assert (dofs != 0, typename DataOut_DoFData<dim>::ExcNoDoFHandlerSelected());
+
+ typedef DataOut_DoFData<dim,dim> BaseClass;
+ Assert (dofs != 0, typename BaseClass::ExcNoDoFHandlerSelected());
#ifdef DEAL_II_USE_MT
const unsigned int n_threads = n_threads_;
// explicit instantiations
-template class DataOut_DoFData<deal_II_dimension>;
+template class DataOut_DoFData<deal_II_dimension,deal_II_dimension>;
+template class DataOut_DoFData<deal_II_dimension,deal_II_dimension+1>;
template class DataOut<deal_II_dimension>;
--- /dev/null
+//---------------------------- data_out_rotation.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_rotation.cc ---------------------------
+
+
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <numerics/data_out_rotation.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>
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void DataOutRotation<1>::build_some_patches (Data)
+{
+ // this function could certainly be
+ // implemented quite easily, but I
+ // haven't done so yet. it may even
+ // be included into the general
+ // template below, if one were not
+ // to choose y as rotational
+ // variable but z, since then one
+ // could loop over all quadrature
+ // points as outer loop and over
+ // the symmetry variable inside, as
+ // the latter is only repeated
+ Assert (false, ExcNotImplemented());
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+void DataOutRotation<2>::build_some_patches (Data data)
+{
+ const unsigned int dim = 2;
+
+ QTrapez<1> q_trapez;
+ QIterated<dim> patch_points (q_trapez, data.n_subdivisions);
+
+ FEValues<dim> fe_patch_values(dofs->get_fe(),
+ patch_points,
+ update_values);
+
+ const unsigned int n_patches_per_circle = data.n_patches_per_circle;
+
+ // another abbreviation denoting
+ // the number of q_points in each
+ // direction
+ const unsigned int n_points = data.n_subdivisions+1;
+
+ // set up an array that holds the
+ // directions in the plane of
+ // rotation in which we will put
+ // points in the whole domain (not
+ // the rotationally reduced one in
+ // which the computation took
+ // place. for simplicity add the
+ // initial direction at the end
+ // again
+ const double pi = 3.14159265358979323846;
+ vector<Point<3> > angle_directions (n_patches_per_circle+1);
+ for (unsigned int i=0; i<=n_patches_per_circle; ++i)
+ {
+ angle_directions[i][0] = cos(2*pi*i/n_patches_per_circle);
+ angle_directions[i][1] = sin(2*pi*i/n_patches_per_circle);
+ };
+
+
+ unsigned int cell_number = 0;
+ vector<DataOutBase::Patch<dim+1> >::iterator patch = patches.begin();
+ DoFHandler<dim>::cell_iterator cell=first_cell();
+
+ // get first cell in this thread
+ for (unsigned int i=0; (i<data.this_thread)&&(cell != dofs->end()); ++i)
+ {
+ advance (patch, n_patches_per_circle);
+ ++cell_number;
+ cell=next_cell(cell);
+ }
+
+ // now loop over all cells and
+ // actually create the patches
+ for (; cell != dofs->end(); )
+ {
+ for (unsigned int angle=0; angle<n_patches_per_circle; ++angle, ++patch)
+ {
+ Assert (patch != patches.end(), ExcInternalError());
+
+
+ // first compute the
+ // vertices of the
+ // patch. note that they
+ // will have to be computed
+ // from the vertices of the
+ // cell, which has one
+ // dimension less, however.
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
+ {
+ const Point<dim> v = cell->vertex(vertex);
+ patch->vertices[vertex] = v(0) * angle_directions[angle];
+ patch->vertices[vertex][2] = v(1);
+
+ patch->vertices[vertex+GeometryInfo<dim>::vertices_per_cell]
+ = v(0) * angle_directions[angle+1];
+ patch->vertices[vertex+GeometryInfo<dim>::vertices_per_cell][2]
+ = v(1);
+ };
+
+ // then fill in data
+ if (data.n_datasets > 0)
+ {
+ fe_patch_values.reinit (cell);
+
+ // 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 x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ for (unsigned int z=0; z<n_points; ++z)
+ patch->data(dataset,
+ x*n_points*n_points +
+ y*n_points +
+ z)
+ = data.patch_values[x*n_points+z];
+ }
+ 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 x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ for (unsigned int z=0; z<n_points; ++z)
+ patch->data(dataset*data.n_components+component,
+ x*n_points*n_points +
+ y*n_points +
+ z)
+ = data.patch_values_system[x*n_points+z](component);
+ };
+ };
+
+ // then do the cell data
+ for (unsigned int dataset=0; dataset<cell_data.size(); ++dataset)
+ {
+ const double value = (*cell_data[dataset].data)(cell_number);
+ for (unsigned int x=0; x<n_points; ++x)
+ for (unsigned int y=0; y<n_points; ++y)
+ for (unsigned int z=0; z<n_points; ++z)
+ patch->data(dataset+dof_data.size()*data.n_components,
+ x*n_points*n_points +
+ y*n_points +
+ z)
+ = value;
+ };
+ };
+ };
+
+ // next cell (patch) in this
+ // thread. note that we have
+ // already advanced the patches
+ // for the present cell,
+ // i.e. we only have to skip
+ // the cells belonging to other
+ // threads, not the ones
+ // belonging to this thread.
+ const int skip_threads = static_cast<signed int>(data.n_threads)-1;
+ for (int i=0; (i<skip_threads) && (cell != dofs->end()); ++i)
+ advance (patch, n_patches_per_circle);
+
+ // however, cell and cell
+ // number have not yet been
+ // increased
+ for (unsigned int i=0; (i<data.n_threads) && (cell != dofs->end()); ++i)
+ {
+ ++cell_number;
+ cell=next_cell(cell);
+ };
+ };
+};
+
+#endif
+
+
+#if deal_II_dimension == 3
+
+template <>
+void DataOutRotation<2>::build_some_patches (Data)
+{
+ // would this function make any
+ // sense after all? who would want
+ // to output/compute in four space
+ // dimensions?
+ Assert (false, ExcNotImplemented());
+};
+
+#endif
+
+
+
+template <int dim>
+void DataOutRotation<dim>::build_patches (const unsigned int n_patches_per_circle,
+ 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> 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> > 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 (DoFHandler<dim>::cell_iterator cell=first_cell();
+ cell != dofs->end();
+ cell = next_cell(cell))
+ ++n_patches;
+ // then also take into account that
+ // we want more than one patch to
+ // come out of every cell, as they
+ // are repeated around the axis of
+ // rotation
+ n_patches *= n_patches_per_circle;
+
+ 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_patches_per_circle = n_patches_per_circle;
+ 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 cell have to be
+ // repeated in angular direction
+ DataOutBase::Patch<dim+1> default_patch;
+ default_patch.n_subdivisions = n_subdivisions;
+ default_patch.data.reinit (n_datasets,
+ n_q_points * (n_subdivisions+1));
+ 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 (&DataOutRotation<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 DoFHandler<dim>::cell_iterator
+DataOutRotation<dim>::first_cell ()
+{
+ return dofs->begin_active ();
+};
+
+
+template <int dim>
+typename DoFHandler<dim>::cell_iterator
+DataOutRotation<dim>::next_cell (const typename DoFHandler<dim>::cell_iterator &cell)
+{
+ // convert the iterator to an
+ // active_iterator and advance
+ // this to the next active cell
+ typename DoFHandler<dim>::active_cell_iterator active_cell = cell;
+ ++active_cell;
+ return active_cell;
+};
+
+
+// explicit instantiations
+template class DataOutRotation<deal_II_dimension>;
+
+
<h3>deal.II</h3>
<ol>
+ <li> <p>
+ New: there is now a class <code class="class">DataOutRotation</code>
+ that can be used to output data which has been computed
+ exploiting rotational symmetry, on the original domain. Thus,
+ the output is of one dimension higher than the computation was,
+ where the computed solution is rotated around the axis of
+ symmetry.
+ <br>
+ (WB 2000/08/14)
+ </p>
+
+ <li> <p>
+ New: class <code class="class">HalfHyperShellBoundary</code>
+ and <code class="member">GridGenerator::half_hyper_shell</code>
+ generate a half shell, useful for computations with a shell
+ domain and rotational symmetry.
+ <br>
+ (WB 2000/08/08)
+ </p>
+
<li> <p>
Changed: The functions
<code class="member">Triangulation::refine</code>,
projects / dealii-svn.git / commitdiff