From: Wolfgang Bangerth Date: Thu, 7 Sep 2000 16:23:51 +0000 (+0000) Subject: Implement a way to output faces instead of cells. X-Git-Tag: v8.0.0~20124 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e0748dc496e07aa1b88c47158ef68f263848e07f;p=dealii.git Implement a way to output faces instead of cells. git-svn-id: https://svn.dealii.org/trunk@3305 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/data_out.h b/deal.II/deal.II/include/numerics/data_out.h index e6349405d0..fef249b50a 100644 --- a/deal.II/deal.II/include/numerics/data_out.h +++ b/deal.II/deal.II/include/numerics/data_out.h @@ -133,8 +133,8 @@ template class DoFHandler; * * @author Wolfgang Bangerth, 1999 */ -template -class DataOut_DoFData : public DataOutInterface +template +class DataOut_DoFData : public DataOutInterface { public: /** @@ -377,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 @@ -385,7 +385,8 @@ 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 @@ -397,6 +398,7 @@ class DataOut_DoFData : public DataOutInterface }; + /** * 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 index 0000000000..16126c4588 --- /dev/null +++ b/deal.II/deal.II/include/numerics/data_out_faces.h @@ -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 + +#include +#include + +template 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 +class DataOutFaces : 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_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{::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::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 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 053bd3a70c..6db618835d 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -28,8 +28,8 @@ #include -template -DataOut_DoFData::DataEntry:: +template +DataOut_DoFData::DataEntry:: DataEntry (const Vector *data, const vector &names) : data(data), @@ -38,24 +38,24 @@ DataEntry (const Vector *data, -template -DataOut_DoFData::DataOut_DoFData () : +template +DataOut_DoFData::DataOut_DoFData () : dofs(0) {}; -template -DataOut_DoFData::~DataOut_DoFData () +template +DataOut_DoFData::~DataOut_DoFData () { clear (); }; -template +template void -DataOut_DoFData:: +DataOut_DoFData:: attach_dof_handler (const DoFHandler &d) { Assert (dof_data.size() == 0, ExcOldDataStillPresent()); @@ -71,9 +71,9 @@ attach_dof_handler (const DoFHandler &d) -template +template void -DataOut_DoFData::add_data_vector (const Vector &vec, +DataOut_DoFData::add_data_vector (const Vector &vec, const vector &names) { Assert (dofs != 0, ExcNoDoFHandlerSelected ()); @@ -112,9 +112,9 @@ DataOut_DoFData::add_data_vector (const Vector +template void -DataOut_DoFData::add_data_vector (const Vector &vec, +DataOut_DoFData::add_data_vector (const Vector &vec, const string &name) { unsigned int n_components = dofs->get_fe().n_components (); @@ -146,22 +146,22 @@ DataOut_DoFData::add_data_vector (const Vector -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 +template void -DataOut_DoFData::clear_input_data_references () +DataOut_DoFData::clear_input_data_references () { for (unsigned int i=0; i::clear_input_data_references () -template +template void -DataOut_DoFData::clear () +DataOut_DoFData::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::clear () }; // delete patches - vector > dummy; + vector > dummy; patches.swap (dummy); } -template +template vector -DataOut_DoFData::get_dataset_names () const +DataOut_DoFData::get_dataset_names () const { vector names; // collect the names of dof @@ -220,9 +220,9 @@ DataOut_DoFData::get_dataset_names () const -template -const vector > & -DataOut_DoFData::get_patches () const +template +const vector > & +DataOut_DoFData::get_patches () const { return patches; }; @@ -430,3 +430,7 @@ template class DataOut_DoFData; template class DataOut; +// for 3d, also generate an extra class +#if deal_II_dimension >= 3 +template class DataOut_DoFData; +#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 index 0000000000..c04b828c60 --- /dev/null +++ b/deal.II/deal.II/source/numerics/data_out_faces.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_MT +#include +#endif + +#include + + + +template +void DataOutFaces::build_some_patches (Data data) +{ + QTrapez<1> q_trapez; + QIterated patch_points (q_trapez, data.n_subdivisions); + + FEFaceValues 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 >::iterator patch = patches.begin(); + FaceDescriptor face=first_face(); + + // get first face in this thread + for (unsigned int i=0; (iend()); ++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::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; datasetdata(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; componentdata(dataset*data.n_components+component,q) = + data.patch_values_system[q](component); + }; + }; + + // then do the cell data + for (unsigned int dataset=0; datasetdata(dataset+dof_data.size()*data.n_components,q) = + value; + }; + }; + // next cell (patch) in this thread + for (unsigned int i=0; + (iend()); ++i) + { + ++patch; + ++face_number; + face=next_face(face); + } + }; +}; + + + + +template +void DataOutFaces::build_patches (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 (FaceDescriptor face=first_face(); + face.first != dofs->end(); + face = next_face(face)) + ++n_patches; + + 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); + 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 DataOutFaces::FaceDescriptor +DataOutFaces::first_face () +{ + // simply find first active cell + // with a face on the boundary + typename DoFHandler::active_cell_iterator cell = dofs->begin_active(); + for (; cell != dofs->end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + return make_pair (static_cast(cell), f); + + // ups, triangulation has no + // boundary? impossible! + Assert (false, ExcInternalError()); + + return FaceDescriptor(); +}; + + + +template +typename DataOutFaces::FaceDescriptor +DataOutFaces::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::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::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::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; +#endif +