From 5369ea155863d859535674605d232d9d690b2ff7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 14 Aug 2000 08:36:37 +0000 Subject: [PATCH] Implement DataOutRotation. git-svn-id: https://svn.dealii.org/trunk@3253 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/data_out.h | 27 +- .../include/numerics/data_out_rotation.h | 128 +++++++ deal.II/deal.II/source/numerics/data_out.cc | 77 ++-- .../source/numerics/data_out_rotation.cc | 356 ++++++++++++++++++ deal.II/doc/news/2000/c-3-0.html | 20 + 5 files changed, 569 insertions(+), 39 deletions(-) create mode 100644 deal.II/deal.II/include/numerics/data_out_rotation.h create mode 100644 deal.II/deal.II/source/numerics/data_out_rotation.cc diff --git a/deal.II/deal.II/include/numerics/data_out.h b/deal.II/deal.II/include/numerics/data_out.h index 9934b77827..e1ecc36634 100644 --- a/deal.II/deal.II/include/numerics/data_out.h +++ b/deal.II/deal.II/include/numerics/data_out.h @@ -118,10 +118,23 @@ template class DoFHandler; * 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 -class DataOut_DoFData : public DataOutInterface +template +class DataOut_DoFData : public DataOutInterface { public: /** @@ -140,7 +153,7 @@ class DataOut_DoFData : public DataOutInterface * and the mapping between nodes * and node values. */ - void attach_dof_handler (const DoFHandler &); + void attach_dof_handler (const DoFHandler &); /** * Add a data vector together @@ -343,7 +356,7 @@ class DataOut_DoFData : public DataOutInterface /** * Pointer to the dof handler object. */ - const DoFHandler *dofs; + const DoFHandler *dofs; /** * List of data elements with vectors of @@ -364,7 +377,7 @@ class DataOut_DoFData : public DataOutInterface * in the output routines of the base * classes. */ - vector > patches; + vector > patches; /** * Function by which the base @@ -372,7 +385,7 @@ class DataOut_DoFData : public DataOutInterface * what patches they shall write * to a file. */ - virtual const vector > & get_patches () const; + virtual const vector > & get_patches () const; /** * Virtual function through @@ -443,7 +456,7 @@ class DataOut_DoFData : public DataOutInterface * @author Wolfgang Bangerth, 1999 */ template -class DataOut : public DataOut_DoFData +class DataOut : public DataOut_DoFData { public: /** diff --git a/deal.II/deal.II/include/numerics/data_out_rotation.h b/deal.II/deal.II/include/numerics/data_out_rotation.h new file mode 100644 index 0000000000..e49e215b6e --- /dev/null +++ b/deal.II/deal.II/include/numerics/data_out_rotation.h @@ -0,0 +1,128 @@ +//---------------------------- 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 + +#include +#include + +template class DoFHandler; + +/** + * + * @author Wolfgang Bangerth, 2000 + */ +template +class DataOutRotation : public DataOut_DoFData +{ + 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::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::cell_iterator + next_cell (const typename DoFHandler::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 patch_values; + vector > 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 diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index 3a5bdb09b9..053bd3a70c 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -28,32 +28,35 @@ #include -template -DataOut_DoFData::DataEntry::DataEntry (const Vector *data, - const vector &names) : +template +DataOut_DoFData::DataEntry:: +DataEntry (const Vector *data, + const vector &names) : data(data), names(names) {}; -template -DataOut_DoFData::DataOut_DoFData () : +template +DataOut_DoFData::DataOut_DoFData () : dofs(0) {}; -template -DataOut_DoFData::~DataOut_DoFData () +template +DataOut_DoFData::~DataOut_DoFData () { clear (); }; -template -void DataOut_DoFData::attach_dof_handler (const DoFHandler &d) +template +void +DataOut_DoFData:: +attach_dof_handler (const DoFHandler &d) { Assert (dof_data.size() == 0, ExcOldDataStillPresent()); Assert (cell_data.size() == 0, ExcOldDataStillPresent()); @@ -68,9 +71,10 @@ void DataOut_DoFData::attach_dof_handler (const DoFHandler &d) -template -void DataOut_DoFData::add_data_vector (const Vector &vec, - const vector &names) +template +void +DataOut_DoFData::add_data_vector (const Vector &vec, + const vector &names) { Assert (dofs != 0, ExcNoDoFHandlerSelected ()); @@ -108,9 +112,10 @@ void DataOut_DoFData::add_data_vector (const Vector &vec, }; -template -void DataOut_DoFData::add_data_vector (const Vector &vec, - const string &name) +template +void +DataOut_DoFData::add_data_vector (const Vector &vec, + const string &name) { unsigned int n_components = dofs->get_fe().n_components (); @@ -141,21 +146,22 @@ void DataOut_DoFData::add_data_vector (const Vector &vec, -template -void DataOut_DoFData::clear_data_vectors () +template +void DataOut_DoFData::clear_data_vectors () { dof_data.erase (dof_data.begin(), dof_data.end()); cell_data.erase (cell_data.begin(), cell_data.end()); // delete patches - vector > dummy; + vector > dummy; patches.swap (dummy); } -template -void DataOut_DoFData::clear_input_data_references () +template +void +DataOut_DoFData::clear_input_data_references () { for (unsigned int i=0; i::clear_input_data_references () - -template -void DataOut_DoFData::clear () +template +void +DataOut_DoFData::clear () { dof_data.erase (dof_data.begin(), dof_data.end()); cell_data.erase (cell_data.begin(), cell_data.end()); @@ -186,13 +192,14 @@ void DataOut_DoFData::clear () }; // delete patches - vector > dummy; + vector > dummy; patches.swap (dummy); } -template -vector DataOut_DoFData::get_dataset_names () const +template +vector +DataOut_DoFData::get_dataset_names () const { vector names; // collect the names of dof @@ -213,15 +220,19 @@ vector DataOut_DoFData::get_dataset_names () const -template -const vector > & -DataOut_DoFData::get_patches () const +template +const vector > & +DataOut_DoFData::get_patches () const { return patches; }; +/* ---------------------------------------------------------------------- */ + + + template void DataOut::build_some_patches (Data data) { @@ -309,8 +320,9 @@ void DataOut::build_patches (const unsigned int n_subdivisions, { Assert (n_subdivisions >= 1, ExcInvalidNumberOfSubdivisions(n_subdivisions)); - - Assert (dofs != 0, typename DataOut_DoFData::ExcNoDoFHandlerSelected()); + + typedef DataOut_DoFData BaseClass; + Assert (dofs != 0, typename BaseClass::ExcNoDoFHandlerSelected()); #ifdef DEAL_II_USE_MT const unsigned int n_threads = n_threads_; @@ -413,7 +425,8 @@ DataOut::next_cell (const typename DoFHandler::cell_iterator &cell) // explicit instantiations -template class DataOut_DoFData; +template class DataOut_DoFData; +template class DataOut_DoFData; template class DataOut; diff --git a/deal.II/deal.II/source/numerics/data_out_rotation.cc b/deal.II/deal.II/source/numerics/data_out_rotation.cc new file mode 100644 index 0000000000..5af090af77 --- /dev/null +++ b/deal.II/deal.II/source/numerics/data_out_rotation.cc @@ -0,0 +1,356 @@ +//---------------------------- 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 +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_MT +#include +#endif + +#include + + + +#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 patch_points (q_trapez, data.n_subdivisions); + + FEValues 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 > 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 >::iterator patch = patches.begin(); + DoFHandler::cell_iterator cell=first_cell(); + + // get first cell in this thread + for (unsigned int i=0; (iend()); ++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::vertices_per_cell; ++vertex) + { + const Point v = cell->vertex(vertex); + patch->vertices[vertex] = v(0) * angle_directions[angle]; + patch->vertices[vertex][2] = v(1); + + patch->vertices[vertex+GeometryInfo::vertices_per_cell] + = v(0) * angle_directions[angle+1]; + patch->vertices[vertex+GeometryInfo::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; datasetdata(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; componentdata(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; datasetdata(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(data.n_threads)-1; + for (int i=0; (iend()); ++i) + advance (patch, n_patches_per_circle); + + // however, cell and cell + // number have not yet been + // increased + for (unsigned int i=0; (iend()); ++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 +void DataOutRotation::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 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 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 > 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::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 thread_data(n_threads); + + // init data for the threads + for (unsigned int i=0;i 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::build_some_patches) + .collect_args (this, thread_data[l])); + thread_manager.wait(); + + // just one thread +#else + build_some_patches(thread_data[0]); +#endif +}; + + +template +typename DoFHandler::cell_iterator +DataOutRotation::first_cell () +{ + return dofs->begin_active (); +}; + + +template +typename DoFHandler::cell_iterator +DataOutRotation::next_cell (const typename DoFHandler::cell_iterator &cell) +{ + // convert the iterator to an + // active_iterator and advance + // this to the next active cell + typename DoFHandler::active_cell_iterator active_cell = cell; + ++active_cell; + return active_cell; +}; + + +// explicit instantiations +template class DataOutRotation; + + diff --git a/deal.II/doc/news/2000/c-3-0.html b/deal.II/doc/news/2000/c-3-0.html index 6833cf6bc9..61acffd2c6 100644 --- a/deal.II/doc/news/2000/c-3-0.html +++ b/deal.II/doc/news/2000/c-3-0.html @@ -391,6 +391,26 @@

deal.II

    +

  1. + New: there is now a class DataOutRotation + 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. +
    + (WB 2000/08/14) +

    + +

  2. + New: class HalfHyperShellBoundary + and GridGenerator::half_hyper_shell + generate a half shell, useful for computations with a shell + domain and rotational symmetry. +
    + (WB 2000/08/08) +

    +

  3. Changed: The functions Triangulation::refine, -- 2.39.5