]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Output of curved patches for DataOutFaces if a Mapping is provided to the build_patch...
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 5 Jun 2009 04:18:11 +0000 (04:18 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 5 Jun 2009 04:18:11 +0000 (04:18 +0000)
git-svn-id: https://svn.dealii.org/trunk@18909 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/data_out_faces.h
deal.II/deal.II/source/numerics/data_out_faces.cc

index c889a67ed7c0540e1a373ac75922d2819c060150..a4cdc00d91ed57dffee3c305128faa790bbbef4d 100644 (file)
@@ -45,10 +45,12 @@ namespace internal
                      const unsigned int n_datasets,
                      const unsigned int n_subdivisions,
                      const std::vector<unsigned int> &n_postprocessor_outputs,
+                     const Mapping<dim,spacedim> &mapping,
                      const FE &finite_elements,
                      const UpdateFlags update_flags);
 
        const dealii::hp::QCollection<dim-1> q_collection;
+       const dealii::hp::MappingCollection<dim,spacedim> mapping_collection;
        dealii::hp::FEFaceValues<dim> x_fe_values;
 
        std::vector<Point<dim> > patch_normals;
@@ -146,6 +148,50 @@ class DataOutFaces : public DataOut_DoFData<DH,DH::dimension-1,
     virtual void
     build_patches (const unsigned int n_subdivisions = 0);
 
+                                    /**
+                                     * Same as above, except that the
+                                     * additional first parameter
+                                     * defines a mapping that is to
+                                     * be used in the generation of
+                                     * output. If
+                                     * <tt>n_subdivisions>1</tt>, the
+                                     * points interior of subdivided
+                                     * patches which originate from
+                                     * cells at the boundary of the
+                                     * domain can be computed using the
+                                     * mapping, i.e. a higher order
+                                     * mapping leads to a
+                                     * representation of a curved
+                                     * boundary by using more
+                                     * subdivisions.
+                                     *
+                                     * Even for non-curved cells the
+                                     * mapping argument can be used
+                                     * for the Eulerian mappings (see
+                                     * class MappingQ1Eulerian) where
+                                     * a mapping is used not only to
+                                     * determine the position of
+                                     * points interior to a cell, but
+                                     * also of the vertices.  It
+                                     * offers an opportunity to watch
+                                     * the solution on a deformed
+                                     * triangulation on which the
+                                     * computation was actually
+                                     * carried out, even if the mesh
+                                     * is internally stored in its
+                                     * undeformed configuration and
+                                     * the deformation is only
+                                     * tracked by an additional
+                                     * vector that holds the
+                                     * deformation of each vertex.
+                                     *
+                                     * @todo The @p mapping argument should be
+                                     * replaced by a hp::MappingCollection in
+                                     * case of a hp::DoFHandler.
+                                     */
+    virtual void build_patches (const Mapping<DH::dimension> &mapping,
+                               const unsigned int n_subdivisions = 0);
+
                                     /**
                                      * Declare a way to describe a
                                      * face which we would like to
index 399a7af9b8917468180ecfb7349b98e99e672074..9ae1ef49f566b0e8ec197dac5e9ea03416c13551 100644 (file)
@@ -41,6 +41,7 @@ namespace internal
                  const unsigned int n_datasets,
                  const unsigned int n_subdivisions,
                  const std::vector<unsigned int> &n_postprocessor_outputs,
+                 const Mapping<dim,spacedim> &mapping,
                  const FE &finite_elements,
                  const UpdateFlags update_flags)
                  :
@@ -52,7 +53,9 @@ namespace internal
                                                  n_postprocessor_outputs,
                                                  finite_elements),
                  q_collection (quadrature),
-                 x_fe_values (this->fe_collection,
+                 mapping_collection (mapping),     
+                 x_fe_values (this->mapping_collection,
+                              this->fe_collection,
                               q_collection,
                               update_flags)
     {}
@@ -84,8 +87,24 @@ build_one_patch (const FaceDescriptor *cell_and_face,
                 internal::DataOutFaces::ParallelData<DH::dimension, DH::dimension> &data,
                 DataOutBase::Patch<DH::dimension-1,DH::space_dimension>  &patch)
 {
+                                  // we use the mapping to transform the
+                                  // vertices. However, the mapping works on
+                                  // cells, not faces, so transform the face
+                                  // vertex to a cell vertex, that to a unit
+                                  // cell vertex and then, finally, that to
+                                  // the mapped vertex. In most cases this
+                                  // complicated procedure will be the
+                                  // identity.
   for (unsigned int vertex=0; vertex<GeometryInfo<DH::dimension-1>::vertices_per_cell; ++vertex)
-    patch.vertices[vertex] = cell_and_face->first->face(cell_and_face->second)->vertex(vertex);
+    patch.vertices[vertex] = data.mapping_collection[0].transform_unit_to_real_cell
+                            (cell_and_face->first,
+                             GeometryInfo<DH::dimension>::unit_cell_vertex
+                             (GeometryInfo<dim>::face_to_cell_vertices
+                              (cell_and_face->second,
+                               vertex,
+                               cell_and_face->first->face_orientation(cell_and_face->second),
+                               cell_and_face->first->face_flip(cell_and_face->second),
+                               cell_and_face->first->face_rotation(cell_and_face->second))));
       
   if (data.n_datasets > 0)
     {
@@ -94,8 +113,25 @@ build_one_patch (const FaceDescriptor *cell_and_face,
        = data.x_fe_values.get_present_fe_values ();
 
       const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
-  
-         
+      
+                                      // store the intermediate points
+      Assert(patch.space_dim==DH::dimension, ExcInternalError());
+      const std::vector<Point<DH::dimension> > & 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)+DH::dimension,
+                       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<DH::dimension; ++i)
+       for (unsigned int q=0; q<n_q_points; ++q)
+         patch.data(patch.data.size(0)-DH::dimension+i,q)=q_points[q][i];
+      
                                       // counter for data records
       unsigned int offset=0;
          
@@ -230,6 +266,15 @@ build_one_patch (const FaceDescriptor *cell_and_face,
 
 template <int dim, class DH>
 void DataOutFaces<dim,DH>::build_patches (const unsigned int n_subdivisions_)
+{
+  build_patches (StaticMappingQ1<DH::dimension>::mapping, n_subdivisions_);
+}
+
+
+
+template <int dim, class DH>
+void DataOutFaces<dim,DH>::build_patches (const Mapping<DH::dimension> &mapping,
+                                         const unsigned int n_subdivisions_)
 {
                                   // Check consistency of redundant
                                   // template parameter
@@ -286,11 +331,14 @@ void DataOutFaces<dim,DH>::build_patches (const unsigned int n_subdivisions_)
   for (unsigned int i=0; i<this->dof_data.size(); ++i)
     if (this->dof_data[i]->postprocessor)
       update_flags |= this->dof_data[i]->postprocessor->get_needed_update_flags();
+  update_flags |= update_quadrature_points;
   
   internal::DataOutFaces::ParallelData<DH::dimension, DH::dimension>
     thread_data (patch_points, n_components, n_datasets,
                 n_subdivisions, 
-                n_postprocessor_outputs, this->dofs->get_fe(),
+                n_postprocessor_outputs,
+                mapping,
+                this->dofs->get_fe(),
                 update_flags);
   DataOutBase::Patch<DH::dimension-1,DH::space_dimension> sample_patch;
   sample_patch.n_subdivisions = n_subdivisions;
@@ -388,7 +436,7 @@ DataOutFaces<dim,DH>::next_face (const FaceDescriptor &old_face)
 
 
 // explicit instantiations
-// don't instantiate anything for the 1d and 2d cases
+// don't instantiate anything for the 1d
 #if deal_II_dimension >=2
 template class DataOutFaces<deal_II_dimension, DoFHandler<deal_II_dimension> >;
 template class DataOutFaces<deal_II_dimension, hp::DoFHandler<deal_II_dimension> >;

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.