From: leicht Date: Tue, 30 Jan 2007 06:57:25 +0000 (+0000) Subject: Extend DataOutBase and DataOut to store mapped (quadrature) points in each patch... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=67135ff95873550d72f5a7f2fb6bb9443e044173;p=dealii-svn.git Extend DataOutBase and DataOut to store mapped (quadrature) points in each patch at the boundary. Use a mapping to obtain function values and thus fix a bug for FE_RaviartThomas and FE_ABF elements, which need such a mapping for the function values. git-svn-id: https://svn.dealii.org/trunk@14392 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/data_out_base.h b/deal.II/base/include/base/data_out_base.h index a801136c6d..dd29714272 100644 --- a/deal.II/base/include/base/data_out_base.h +++ b/deal.II/base/include/base/data_out_base.h @@ -108,6 +108,10 @@ class ParameterHandler; * n_subdivisions==3 will yield 4 times 4 (times 4) points, etc. The actual * location of these points on the patch will be computed by a multilinear * transformation from the vertices given for this patch. + + * For cells at the boundary, a mapping might be used to calculate the position + * of the inner points. In that case the coordinates are stored inside the + * Patch, as they cannot be easily recovered otherwise. * * Given these comments, the actual data to be printed on this patch of * points consists of several data sets each of which has a value at each @@ -252,6 +256,12 @@ class DataOutBase template struct Patch { + /** + * Make the spacedim template + * parameter available. + */ + static const unsigned int space_dim=spacedim; + /** * Corner points of a patch. * Inner points are computed by @@ -329,6 +339,23 @@ class DataOutBase * patches provided. */ Table<2,float> data; + + /** + * Bool flag indicating, whether the + * coordinates of the inner patch + * points are appended to the @p data + * table (@ true) or not (@ false), + * where the second is the standard and + * can be found for all cells in the + * interior of a domain. On the + * boundary, patch points are evaluated + * using a Mapping and therefore have + * to be stored inside the patch, as + * the Mapping and the corresponding + * boundary information might not be + * available later on. + */ + bool points_are_available; /** * Default constructor. Sets diff --git a/deal.II/base/source/data_out_base.cc b/deal.II/base/source/data_out_base.cc index 194cd8a1b2..4ead59c4c1 100644 --- a/deal.II/base/source/data_out_base.cc +++ b/deal.II/base/source/data_out_base.cc @@ -96,29 +96,60 @@ static unsigned int vtk_cell_type[4] = //Auxiliary functions //----------------------------------------------------------------------// //For a given patch, compute the node interpolating the corner nodes -//linearly at the point (xfrac, yfrac, zfrac). +//linearly at the point (xstep, ystep, zstep)*1./n_subdivisions. +//If the popints are saved in the patch->data member, return the +//saved point instead template inline void compute_node( Point& node, const DataOutBase::Patch* patch, - const float xfrac, - const float yfrac, - const float zfrac) + const unsigned int xstep, + const unsigned int ystep, + const unsigned int zstep, + const unsigned int n_subdivisions) { - node = (patch->vertices[1] * xfrac) + (patch->vertices[0] * (1-xfrac)); - if (dim>1) + if (patch->points_are_available) { - node*= 1-yfrac; - node += ((patch->vertices[3] * xfrac) + (patch->vertices[2] * (1-xfrac))) * yfrac; - if (dim>2) + unsigned int point_no=0; + // note: switch without break ! + switch (dim) { - node *= (1-zfrac); - node += (((patch->vertices[5] * xfrac) + (patch->vertices[4] * (1-xfrac))) - * (1-yfrac) + - ((patch->vertices[7] * xfrac) + (patch->vertices[6] * (1-xfrac))) - * yfrac) * zfrac; + case 3: + Assert (zstepdata(patch->data.size(0)-spacedim+d,point_no); + } + else + { + // perform a dim-linear interpolation + const double stepsize=1./n_subdivisions, + xfrac=xstep*stepsize; + + node = (patch->vertices[1] * xfrac) + (patch->vertices[0] * (1-xfrac)); + if (dim>1) + { + const double yfrac=ystep*stepsize; + node*= 1-yfrac; + node += ((patch->vertices[3] * xfrac) + (patch->vertices[2] * (1-xfrac))) * yfrac; + if (dim>2) + { + const double zfrac=zstep*stepsize; + node *= (1-zfrac); + node += (((patch->vertices[5] * xfrac) + (patch->vertices[4] * (1-xfrac))) + * (1-yfrac) + + ((patch->vertices[7] * xfrac) + (patch->vertices[6] * (1-xfrac))) + * yfrac) * zfrac; + } } } } @@ -170,6 +201,8 @@ compute_sizes(const std::vector >& patches, //----------------------------------------------------------------------// +template +const unsigned int DataOutBase::Patch::space_dim; @@ -181,7 +214,8 @@ template DataOutBase::Patch::Patch () : patch_index(no_neighbor), - n_subdivisions (1) + n_subdivisions (1), + points_are_available(false) // all the other data has a // constructor of its own, except // for the "neighbors" field, which @@ -215,6 +249,9 @@ DataOutBase::Patch::operator == (const Patch &patch) const if (n_subdivisions != patch.n_subdivisions) return false; + if (points_are_available != patch.points_are_available) + return false; + if (data.n_rows() != patch.data.n_rows()) return false; @@ -240,7 +277,9 @@ DataOutBase::Patch::memory_consumption () const + MemoryConsumption::memory_consumption(n_subdivisions) + - MemoryConsumption::memory_consumption(data)); + MemoryConsumption::memory_consumption(data) + + + MemoryConsumption::memory_consumption(points_are_available)); } @@ -1218,9 +1257,10 @@ DataOutBase::write_nodes ( for (unsigned int i1=0; i1n_subdivisions; const unsigned int n = n_subdivisions+1; // Length of loops in all dimensions - Assert (patch->data.n_rows() == n_data_sets, - ExcDimensionMismatch (patch->data.n_rows(), n_data_sets)); + Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) || + (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available), + ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets)); Assert (patch->data.n_cols() == template_power(n), ExcInvalidDatasetSize (patch->data.n_cols(), n)); @@ -1729,8 +1770,9 @@ void DataOutBase::write_gnuplot (const std::vector > &patche unsigned int d2 = n; unsigned int d3 = n*n; - Assert (patch->data.n_rows() == n_data_sets, - ExcDimensionMismatch (patch->data.n_rows(), n_data_sets)); + Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) || + (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available), + ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets)); Assert (patch->data.n_cols() == template_power(n), ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1)); @@ -1742,12 +1784,9 @@ void DataOutBase::write_gnuplot (const std::vector > &patche { for (unsigned int i1=0; i1 > &patche for (unsigned int i2=0; i2 > &patche // write point there // and its data - const double x_frac_new = x_frac + 1./n_subdivisions; - compute_node(node, &*patch, x_frac_new, y_frac, z_frac); + compute_node(node, &*patch, i1+1, i2, i3, n_subdivisions); out << node; for (unsigned int data_set=0; data_set > &patche // write point there // and its data - const double y_frac_new = y_frac + 1./n_subdivisions; - compute_node(node, &*patch, x_frac, y_frac_new, z_frac); + compute_node(node, &*patch, i1, i2+1, i3, n_subdivisions); out << node; for (unsigned int data_set=0; data_set > &patche // write point there // and its data - const double z_frac_new = z_frac + 1./n_subdivisions; - compute_node(node, &*patch, x_frac, y_frac, z_frac_new); + compute_node(node, &*patch, i1, i2, i3+1, n_subdivisions); out << node; for (unsigned int data_set=0; data_set > &patches { const unsigned int n_subdivisions = patch->n_subdivisions; - Assert (patch->data.n_rows() == n_data_sets, - ExcDimensionMismatch (patch->data.n_rows(), n_data_sets)); + Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) || + (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available), + ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets)); Assert (patch->data.n_cols() == template_power(n_subdivisions+1), ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1)); @@ -2005,8 +2038,9 @@ void DataOutBase::write_povray (const std::vector > &patches const unsigned int d1=1; const unsigned int d2=n; - Assert (patch->data.n_rows() == n_data_sets, - ExcDimensionMismatch (patch->data.n_rows(), n_data_sets)); + Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) || + (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available), + ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets)); Assert (patch->data.n_cols() == template_power(n), ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1)); @@ -2016,11 +2050,9 @@ void DataOutBase::write_povray (const std::vector > &patches for (unsigned int i2=0; i2 > &patches, for (unsigned int i2=0; i2 points[4]; - compute_node(points[0], &*patch, x_frac, y_frac, 0.); - compute_node(points[1], &*patch, x_frac1, y_frac, 0.); - compute_node(points[2], &*patch, x_frac, y_frac1, 0.); - compute_node(points[3], &*patch, x_frac1, y_frac1, 0.); + compute_node(points[0], &*patch, i1, i2, 0, n_subdivisions); + compute_node(points[1], &*patch, i1+1, i2, 0, n_subdivisions); + compute_node(points[2], &*patch, i1, i2+1, 0, n_subdivisions); + compute_node(points[3], &*patch, i1+1, i2+1, 0, n_subdivisions); switch (spacedim) { @@ -2596,8 +2622,9 @@ void DataOutBase::write_gmv (const std::vector > &patches, // first patch. checks against all // other patches are made in // write_gmv_reorder_data_vectors - Assert (n_data_sets == patches[0].data.n_rows(), - ExcDimensionMismatch (patches[0].data.n_rows(), n_data_sets)); + Assert ((patches[0].data.n_rows() == n_data_sets && !patches[0].points_are_available) || + (patches[0].data.n_rows() == n_data_sets+spacedim && patches[0].points_are_available), + ExcDimensionMismatch (patches[0].points_are_available ? (patches[0].data.n_rows() + spacedim) : patches[0].data.n_rows(), n_data_sets)); /////////////////////// @@ -2727,11 +2754,12 @@ void DataOutBase::write_tecplot (const std::vector > &patche // first patch. checks against all // other patches are made in // write_gmv_reorder_data_vectors - Assert (n_data_sets == patches[0].data.n_rows(), - ExcDimensionMismatch (patches[0].data.n_rows(), n_data_sets)); - + Assert ((patches[0].data.n_rows() == n_data_sets && !patches[0].points_are_available) || + (patches[0].data.n_rows() == n_data_sets+spacedim && patches[0].points_are_available), + ExcDimensionMismatch (patches[0].points_are_available ? (patches[0].data.n_rows() + spacedim) : patches[0].data.n_rows(), n_data_sets)); + // first count the number of cells // and cells for later use unsigned int n_nodes; @@ -2994,11 +3022,12 @@ void DataOutBase::write_tecplot_binary (const std::vector > // first patch. checks against all // other patches are made in // write_gmv_reorder_data_vectors - Assert (n_data_sets == patches[0].data.n_rows(), - ExcDimensionMismatch (patches[0].data.n_rows(), n_data_sets)); - + Assert ((patch[0].data.n_rows() == n_data_sets && !patches[0].points_are_available) || + (patches[0].data.n_rows() == n_data_sets+spacedim && patches[0].points_are_available), + ExcDimensionMismatch (patches[0].points_are_available ? (patches[0].data.n_rows() + spacedim) : patches[0].data.n_rows(), n_data_sets)); + // first count the number of cells // and cells for later use unsigned int n_nodes; @@ -3282,8 +3311,9 @@ void DataOutBase::write_vtk (const std::vector > &patches, // first patch. checks against all // other patches are made in // write_gmv_reorder_data_vectors - Assert (n_data_sets == patches[0].data.n_rows(), - ExcDimensionMismatch (patches[0].data.n_rows(), n_data_sets)); + Assert ((patches[0].data.n_rows() == n_data_sets && !patches[0].points_are_available) || + (patches[0].data.n_rows() == n_data_sets+spacedim && patches[0].points_are_available), + ExcDimensionMismatch (patches[0].points_are_available ? (patches[0].data.n_rows() + spacedim) : patches[0].data.n_rows(), n_data_sets)); /////////////////////// @@ -3469,7 +3499,12 @@ DataOutBase::write_gmv_reorder_data_vectors (const std::vectordata table + const unsigned int n_data_sets + =patches[0].points_are_available ? (patches[0].data.n_rows() - spacedim) : patches[0].data.n_rows(); Assert (data_vectors.size()[0] == n_data_sets, ExcInternalError()); @@ -3481,8 +3516,9 @@ DataOutBase::write_gmv_reorder_data_vectors (const std::vectorn_subdivisions; - Assert (patch->data.n_rows() == n_data_sets, - ExcDimensionMismatch (patch->data.n_rows(), n_data_sets)); + Assert ((patch->data.n_rows() == n_data_sets && !patch->points_are_available) || + (patch->data.n_rows() == n_data_sets+spacedim && patch->points_are_available), + ExcDimensionMismatch (patch->points_are_available ? (patch->data.n_rows() + spacedim) : patch->data.n_rows(), n_data_sets)); Assert (patch->data.n_cols() == template_power(n_subdivisions+1), ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1)); @@ -4044,6 +4080,8 @@ operator << (std::ostream &out, out << patch.patch_index << ' ' << patch.n_subdivisions << '\n'; + out << patch.points_are_available<<'\n'; + out << patch.data.n_rows() << ' ' << patch.data.n_cols() << '\n'; for (unsigned int i=0; i> (std::istream &in, for (unsigned int i=0; i::faces_per_cell; ++i) in >> patch.neighbors[i]; - in >> patch.patch_index >> patch.n_subdivisions; + in >> patch.patch_index >> patch.n_subdivisions >> patch.points_are_available; unsigned int n_rows, n_cols; in >> n_rows >> n_cols; diff --git a/deal.II/deal.II/include/numerics/data_out.h b/deal.II/deal.II/include/numerics/data_out.h index 32857b8a0a..ff63d0d4e3 100644 --- a/deal.II/deal.II/include/numerics/data_out.h +++ b/deal.II/deal.II/include/numerics/data_out.h @@ -907,18 +907,18 @@ class DataOut : public DataOut_DoFData * additional first parameter * defines a mapping that is to * be used in the generation of - * output. For internal reasons, - * even if + * output. If * n_subdivisions>1, the * points interior of subdivided - * patches are still computed - * from the vertices of the cell, - * i.e. even a higher order - * mapping does not lead to a + * patches which originate from + * cells at the boundary of the + * domain are computed using the + * mapping, i.e. a higher order + * mapping leads to a * representation of a curved * boundary by using more - * subdivisions. However, the - * mapping argument can be used + * subdivisions. The + * mapping argument can also be used * for the Eulerian mappings (see * class MappingQ1Eulerian) where * a mapping is used not only to diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index a38984342d..25e72cfc1a 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -384,26 +384,28 @@ void DataOut::build_some_patches (Data &data) QTrapez<1> q_trapez; QIterated patch_points (q_trapez, data.n_subdivisions); -//TODO[?]: This is strange -- Data has a member 'mapping' that should -//be used here, but it isn't. Rather, up until version 1.94, we were -//actually initializing a local mapping object and used that... While -//we use the mapping to transform the vertex coordinates, we do not -//use the mapping to transform the shape functions (necessary for -//example for Raviart-Thomas elements). This could lead to trouble -//when someone tries to use MappingEulerian with such elements +// We use the mapping to transform the vertex coordinates and the shape +// functions (necessary for example for Raviart-Thomas elements). On the +// boundary, general mappings do not reduce to a MappingQ1, therefore the mapped +// (quadrature) points are stored in the patch, whereas for cells in the +// interior of the domain these points are obtained by a dim-linear mapping and +// can be recovered from the vertices later on, thus they need not to be stored. // create collection objects from // single quadratures, - // and finite elements. if we have + // mappings and finite elements. if we have // an hp DoFHandler, // dof_handler.get_fe() returns a // collection of which we do a // shallow copy instead const hp::QCollection q_collection (patch_points); const hp::FECollection fe_collection(this->dofs->get_fe()); + const hp::MappingCollection mapping_collection(*(data.mapping)); - hp::FEValues x_fe_patch_values (fe_collection, q_collection, - update_values); + hp::FEValues x_fe_patch_values (mapping_collection, + fe_collection, + q_collection, + update_values | update_q_points); const unsigned int n_q_points = patch_points.n_quadrature_points; @@ -436,6 +438,29 @@ void DataOut::build_some_patches (Data &data) const FEValues &fe_patch_values = x_fe_patch_values.get_present_fe_values (); + // if the cell is at the boundary, + // append the quadrature points to + // the last rows of the + // patch->data member + if (cell->at_boundary()) + { + Assert(patch->space_dim==dim, ExcInternalError()); + const std::vector > & q_points=fe_patch_values.get_quadrature_points(); + // resize the patch->data member + // in order to have enough memory + // for the quadrature points as + // well + patch->data.reinit(patch->data.size(0)+dim,patch->data.size(1)); + // set the flag indicating that + // for this cell the points are + // explicitly given + patch->points_are_available=true; + // copy points to patch->data + for (unsigned int i=0; idata(patch->data.size(0)-dim+i,q)=q_points[q][i]; + } + // first fill dof_data for (unsigned int dataset=0; datasetdof_data.size(); ++dataset) { diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index cce5704b4f..077d3618ab 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -212,6 +212,25 @@ inconvenience this causes.
    +
  1. Extended: DataOutBase::Patch has been extended by a new + boolean flag points_are_available, which defaults to + false. It is set to true if the coordinates of the points + defining the subdivision of a patch are appended to the data table + contained in a Patch. This way, DataOut::build_patches() can + use a Mapping to represent curved boundaries, + especially for higher order elements. This change corresponds to an extension + of the intermediate format for graphics. +
    + Using the given Mapping to obtain function values + in DataOut::build_patches() also fixes a bug for FE_RaviartThomas and FE_ABF + elements, which need to evaluate the function values on the real (mapped) cell. +
    + (Tobias Leicht 2007/01/30) +

    +
  2. New: a program reconfigure has been added to the deal.II main directory, which reruns configure with the sam command line arguments as last time.