]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Extend DataOutBase and DataOut to store mapped (quadrature) points in each patch...
authorleicht <leicht@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Jan 2007 06:57:25 +0000 (06:57 +0000)
committerleicht <leicht@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Jan 2007 06:57:25 +0000 (06:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@14392 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/data_out_base.h
deal.II/base/source/data_out_base.cc
deal.II/deal.II/include/numerics/data_out.h
deal.II/deal.II/source/numerics/data_out.cc
deal.II/doc/news/changes.html

index a801136c6d40b536a7b29dc9c593314239751a65..dd297142725a253129f11b8f549f55d62ed836f9 100644 (file)
@@ -108,6 +108,10 @@ class ParameterHandler;
  * <tt>n_subdivisions==3</tt> 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 <int dim, int spacedim=dim>
     struct Patch
     {
+                                        /**
+                                         * Make the <tt>spacedim</tt> 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
index 194cd8a1b29b7055f9c1d7d988cc6af10563a156..4ead59c4c185aba8b3fcc4aa5f4f9284512f1bdd 100644 (file)
@@ -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 <int dim, int spacedim>
 inline
 void
 compute_node(
   Point<spacedim>& node,
   const DataOutBase::Patch<dim,spacedim>* 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 (zstep<n_subdivisions+1, ExcIndexRange(zstep,0,n_subdivisions+1));
+               point_no+=(n_subdivisions+1)*(n_subdivisions+1)*zstep;
+         case 2:
+               Assert (ystep<n_subdivisions+1, ExcIndexRange(ystep,0,n_subdivisions+1));
+               point_no+=(n_subdivisions+1)*ystep;
+         case 1:
+               Assert (xstep<n_subdivisions+1, ExcIndexRange(xstep,0,n_subdivisions+1));
+               point_no+=xstep;
+       }
+      for (unsigned int d=0; d<spacedim; ++d)
+       node[d]=patch->data(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<DataOutBase::Patch<dim, spacedim> >& patches,
 
   
 //----------------------------------------------------------------------//
+template <int dim, int spacedim>
+const unsigned int DataOutBase::Patch<dim,spacedim>::space_dim;
 
 
 
@@ -181,7 +214,8 @@ template <int dim, int spacedim>
 DataOutBase::Patch<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::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; i1<n1; ++i1)
            {
              compute_node(node, &*patch,
-                          i1 * 1./n_subdivisions,
-                          i2 * 1./n_subdivisions,
-                          i3 * 1./n_subdivisions);
+                          i1,
+                          i2,
+                          i3,
+                          n_subdivisions);
              out.write(count++, node);
            }
     }
@@ -1286,8 +1326,9 @@ DataOutBase::write_data (
       const unsigned int n_subdivisions = patch->n_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<dim>(n),
              ExcInvalidDatasetSize (patch->data.n_cols(), n));
       
@@ -1729,8 +1770,9 @@ void DataOutBase::write_gnuplot (const std::vector<Patch<dim,spacedim> > &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<dim>(n),
              ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1));
 
@@ -1742,12 +1784,9 @@ void DataOutBase::write_gnuplot (const std::vector<Patch<dim,spacedim> > &patche
            {
              for (unsigned int i1=0; i1<n1; ++i1)
                {
-                 const double x_frac = i1 * 1./n_subdivisions,
-                              y_frac = i2 * 1./n_subdivisions;
-                 
                                                   // compute coordinates for
                                                   // this patch point
-                 compute_node(node, &*patch, x_frac, y_frac, 0.);
+                 compute_node(node, &*patch, i1, i2, 0, n_subdivisions);
                  out << node << ' ';
                  
                  for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
@@ -1774,13 +1813,9 @@ void DataOutBase::write_gnuplot (const std::vector<Patch<dim,spacedim> > &patche
            for (unsigned int i2=0; i2<n2; ++i2)
              for (unsigned int i1=0; i1<n1; ++i1)
                {
-                 const double x_frac = i1 * 1./n_subdivisions,
-                              y_frac = i2 * 1./n_subdivisions,
-                              z_frac = i3 * 1./n_subdivisions;
-                 
                                                   // compute coordinates for
                                                   // this patch point
-                 compute_node(this_point, &*patch, x_frac, y_frac, z_frac);
+                 compute_node(this_point, &*patch, i1, i2, i3, n_subdivisions);
                                                   // line into positive x-direction
                                                   // if possible
                  if (i1 < n_subdivisions)
@@ -1795,8 +1830,7 @@ void DataOutBase::write_gnuplot (const std::vector<Patch<dim,spacedim> > &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<n_data_sets; ++data_set)
@@ -1823,8 +1857,7 @@ void DataOutBase::write_gnuplot (const std::vector<Patch<dim,spacedim> > &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<n_data_sets; ++data_set)
@@ -1851,8 +1884,7 @@ void DataOutBase::write_gnuplot (const std::vector<Patch<dim,spacedim> > &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<n_data_sets; ++data_set)
@@ -1956,8 +1988,9 @@ void DataOutBase::write_povray (const std::vector<Patch<dim,spacedim> > &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<dim>(n_subdivisions+1),
              ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1));
       
@@ -2005,8 +2038,9 @@ void DataOutBase::write_povray (const std::vector<Patch<dim,spacedim> > &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<dim>(n),
              ExcInvalidDatasetSize (patch->data.n_cols(), n_subdivisions+1));
 
@@ -2016,11 +2050,9 @@ void DataOutBase::write_povray (const std::vector<Patch<dim,spacedim> > &patches
       for (unsigned int i2=0; i2<n; ++i2)
        for (unsigned int i1=0; i1<n; ++i1)
          {
-           const double x_frac = i1 * 1./n_subdivisions,
-                        y_frac = i2 * 1./n_subdivisions;
                                             // compute coordinates for
                                             // this patch point, storing in ver
-           compute_node(ver[i1*d1+i2*d2], &*patch, x_frac, y_frac, 0.);
+           compute_node(ver[i1*d1+i2*d2], &*patch, i1, i2, 0, n_subdivisions);
          }
 
       
@@ -2271,17 +2303,11 @@ void DataOutBase::write_eps (const std::vector<Patch<dim,spacedim> > &patches,
       for (unsigned int i2=0; i2<n_subdivisions; ++i2)
        for (unsigned int i1=0; i1<n_subdivisions; ++i1)
          {
-           const double x_frac = i1 * 1./n_subdivisions,
-                        y_frac = i2 * 1./n_subdivisions,
-                                 
-                       x_frac1 = (i1+1) * 1./n_subdivisions,
-                       y_frac1 = (i2+1) * 1./n_subdivisions;
-           
            Point<spacedim> 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<Patch<dim,spacedim> > &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<Patch<dim,spacedim> > &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<Patch<dim,spacedim> >
                                   // 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<Patch<dim,spacedim> > &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::vector<Patch<dim,spacedi
                                   // first patch. the equivalence of
                                   // these two definitions is checked
                                   // in the main function.
-  const unsigned int n_data_sets = patches[0].data.n_rows();
+  
+                                  // we have to take care, however, whether the
+                                  // points are appended to the end of the
+                                  // patch->data 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::vector<Patch<dim,spacedi
     {
       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<dim>(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<patch.data.n_rows(); ++i)
     for (unsigned int j=0; j<patch.data.n_cols(); ++j)
@@ -4093,7 +4131,7 @@ operator >> (std::istream                     &in,
   for (unsigned int i=0; i<GeometryInfo<dim>::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;
index 32857b8a0ad3dcc74fb9e886fc8d2b26b4fe02c3..ff63d0d4e3005e0b91f53708893f11ce2cbfc296 100644 (file)
@@ -907,18 +907,18 @@ class DataOut : public DataOut_DoFData<DH, DH::dimension>
                                      * additional first parameter
                                      * defines a mapping that is to
                                      * be used in the generation of
-                                     * output. For internal reasons,
-                                     * even if
+                                     * output. If
                                      * <tt>n_subdivisions>1</tt>, 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
index a38984342d5bf12aa4ba7b8a85f9303833d46787..25e72cfc1a3f719a095f1ca5acbbdbc5b17de70f 100644 (file)
@@ -384,26 +384,28 @@ void DataOut<dim,DH>::build_some_patches (Data &data)
   QTrapez<1>     q_trapez;
   QIterated<DH::dimension> 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<DH::dimension>       q_collection (patch_points);
   const hp::FECollection<DH::dimension>      fe_collection(this->dofs->get_fe());
+  const hp::MappingCollection<DH::dimension> mapping_collection(*(data.mapping));
   
-  hp::FEValues<DH::dimension> x_fe_patch_values (fe_collection, q_collection,
-                                       update_values);
+  hp::FEValues<DH::dimension> 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<dim,DH>::build_some_patches (Data &data)
           const FEValues<DH::dimension> &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<Point<dim> > & 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; i<dim; ++i)
+               for (unsigned int q=0; q<n_q_points; ++q)
+                 patch->data(patch->data.size(0)-dim+i,q)=q_points[q][i];
+           }
+
                                           // first fill dof_data
          for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
            {
index cce5704b4ff455a268e831bb31df3123e18b16d0..077d3618abfcfe58e0b9d6b84040112e457c25fc 100644 (file)
@@ -212,6 +212,25 @@ inconvenience this causes.
 
 <ol>
 
+  <li> <p>Extended: <code>DataOutBase::Patch</code> has been extended by a new
+  boolean flag <tt>points_are_available</tt>, which defaults to
+  <tt>false</tt>. It is set to <tt>true</tt> if the coordinates of the points
+  defining the subdivision of a patch are appended to the <tt>data</tt> table
+  contained in a <code>Patch</code>. This way, <code
+  class="class">DataOut</code>::<code class="member">build_patches()</code> can
+  use a <code class="class">Mapping</code> to represent curved boundaries,
+  especially for higher order elements. This change corresponds to an extension
+  of the intermediate format for graphics.
+  <br>
+  Using the given <code class="class">Mapping</code> to obtain function values
+  in <code class="class">DataOut</code>::<code
+  class="member">build_patches()</code> also fixes a bug for <code
+  class="class">FE_RaviartThomas</code> and <code class="class">FE_ABF</code>
+  elements, which need to evaluate the function values on the real (mapped) cell. 
+  <br>
+  (Tobias Leicht 2007/01/30)
+  </p>
+
   <li> <p>New: a program <tt>reconfigure</tt> has been added to the
   deal.II main directory, which reruns <tt>configure</tt> with the sam
   command line arguments as last time.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.