]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a bug in the MappingQ::transform_unit_to_real(real_to_unit)_cell functions. This...
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 28 Mar 2001 15:03:15 +0000 (15:03 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 28 Mar 2001 15:03:15 +0000 (15:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@4316 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/mapping_q.h
deal.II/deal.II/include/fe/mapping_q1.h
deal.II/deal.II/source/fe/mapping_q.cc
deal.II/deal.II/source/fe/mapping_q1.cc

index 7b897388fa2a063442ebaa0f174c109b2d3e2204..d894a80537dbf9164e491cb77ace5d31a0bb913c 100644 (file)
@@ -50,6 +50,30 @@ class MappingQ : public MappingQ1<dim>
                                      * Destructor.
                                      */
     virtual ~MappingQ ();
+    
+                                    /**
+                                     * Transforms the point @p{p} on
+                                     * the unit cell to the point
+                                     * @p{p_real} on the real cell
+                                     * @p{cell} and returns @p{p_real}.
+                                     */
+    virtual Point<dim> transform_unit_to_real_cell (
+      const typename Triangulation<dim>::cell_iterator cell,
+      const Point<dim> &p) const;
+    
+                                    /**
+                                     * Transforms the point @p{p} on
+                                     * the real cell to the point
+                                     * @p{p_unit} on the unit cell
+                                     * @p{cell} and returns @p{p_unit}.
+                                     *
+                                     * Uses Newton iteration and the
+                                     * @p{transform_unit_to_real_cell}
+                                     * function.
+                                     */
+    virtual Point<dim> transform_real_to_unit_cell (
+      const typename Triangulation<dim>::cell_iterator cell,
+      const Point<dim> &p) const;
 
                                     /**
                                      * Implementation of the interface in
index 49e59156e2af4324d0dcf2324a37d6d8f93a3f24..7d15d8310a4314b7f959c20a0d3464c22114cb42 100644 (file)
@@ -388,15 +388,6 @@ class MappingQ1 : public Mapping<dim>
     virtual void compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
                                         InternalData &data) const;
 
-                                    /**
-                                     * Mapping between tensor product
-                                     * ordering and real ordering of
-                                     * the vertices.
-                                     */
-    static const unsigned int vertex_mapping[1<<dim];
-    
-
-  private:
                                     /**
                                      * Transforms a point @p{p} on
                                      * the unit cell to the point
@@ -407,7 +398,7 @@ class MappingQ1 : public Mapping<dim>
                                      * @p{transform_unit_to_real_cell}
                                      * and multiply (through the
                                      * Newton iteration) by
-                                     * @p{transform_real_to_unit_cell}.
+                                     * @p{transform_real_to_unit_cell_internal}.
                                      *
                                      * Takes a reference to an
                                      * @p{InternalData} that must
@@ -423,31 +414,55 @@ class MappingQ1 : public Mapping<dim>
                                      * computations of the mapping
                                      * support points.
                                      */
-    virtual Point<dim> transform_unit_to_real_cell_internal (const InternalData &m_data) const;
-
+    Point<dim> transform_unit_to_real_cell_internal (const InternalData &mdata) const;
+    
                                     /**
-                                     * Returns an @p{InternalData}
-                                     * whose data vectors are resized
-                                     * corresponding to the
-                                     * @p{update_flags} and a
-                                     * one-point
-                                     * quadrature. Furthermore the
-                                     * @p{InternalData} stores the
-                                     * mapping support points of the
-                                     * given @p{cell}.
+                                     * Transforms the point @p{p} on
+                                     * the real cell to the point
+                                     * @p{p_unit} on the unit cell
+                                     * @p{cell} by a Newton
+                                     * iteration.
                                      *
-                                     * This function is called by
-                                     * @p{transform_unit_to_real_cell}
-                                     * and by
-                                     * @p{transform_real_to_unit_cell}.
-                                     * The resulting @p{InternalData}
-                                     * is given to the
-                                     * @p{transform_unit_to_real_internal}
-                                     * function.
+                                     * Takes a reference to an
+                                     * @p{InternalData} that is
+                                     * assumed to be previously
+                                     * created by the @p{get_data}
+                                     * function with @p{UpdateFlags}
+                                     * including
+                                     * @p{update_transformation_values}
+                                     * and
+                                     * @p{update_transformation_gradients}
+                                     * and a one point Quadrature
+                                     * including the given point
+                                     * @p{p_unit}.  Hence this
+                                     * function assumes that
+                                     * @p{mdata} already includes the
+                                     * transformation shape values
+                                     * and gradients computed at
+                                     * @p{p_unit}.
+                                     *
+                                     * These assumptions should be
+                                     * fulfilled by the calling
+                                     * function. That is up to now
+                                     * only the function
+                                     * @p{transform_real_to_unit_cell}
+                                     * and its overloaded versions.
+                                     * @p{mdata} will be changed by
+                                     * this function.
                                      */
-    InternalData* get_cell_data(const typename Triangulation<dim>::cell_iterator cell,
-                               const UpdateFlags update_flags) const;
+    void transform_real_to_unit_cell_internal (const typename Triangulation<dim>::cell_iterator cell,
+                                              const Point<dim> &p,
+                                              InternalData &mdata,
+                                              Point<dim> &p_unit) const;
     
+                                    /**
+                                     * Mapping between tensor product
+                                     * ordering and real ordering of
+                                     * the vertices.
+                                     */
+    static const unsigned int vertex_mapping[1<<dim];
+
+  private:
                                     /**
                                      * Implementation of the interface in
                                      * @ref{Mapping}.
index 365f26c6ee6d6de87666b24ddd0e6dbdade8bedd..a1004e0741cd7ac3816c825757b14b1689a5dd03 100644 (file)
@@ -147,7 +147,6 @@ void
 MappingQ<dim>::compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
                                       MappingQ1<dim>::InternalData &data) const
 {
-  
   const unsigned int n_points=unit_points.size();
   std::vector<double> values;
   std::vector<Tensor<1,dim> > grads;
@@ -1208,6 +1207,76 @@ MappingQ<dim>::transform_contravariant (std::vector<Point<dim> >       &dst,
 }
 
 
+template <int dim>
+Point<dim> MappingQ<dim>::transform_unit_to_real_cell (
+  const typename Triangulation<dim>::cell_iterator cell,
+  const Point<dim> &p) const
+{
+                                  // Use the get_data function to
+                                  // create an InternalData with data
+                                  // vectors of the right size and
+                                  // transformation shape values
+                                  // already computed at point p.
+  const Quadrature<dim> point_quadrature(p);
+  InternalData *mdata=dynamic_cast<InternalData *> (
+    get_data(update_transformation_values, point_quadrature));
+  Assert(mdata!=0, ExcInternalError());
+  
+  mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
+                                            || cell->has_boundary_lines());
+
+  MappingQ1<dim>::InternalData *p_data=0;
+  if (mdata->use_mapping_q1_on_current_cell)
+    p_data=&mdata->mapping_q1_data;
+  else
+    p_data=mdata;
+
+  compute_mapping_support_points(cell, p_data->mapping_support_points);
+  
+  return transform_unit_to_real_cell_internal(*p_data);
+}
+
 
+template <int dim>
+Point<dim> MappingQ<dim>::transform_real_to_unit_cell (
+  const typename Triangulation<dim>::cell_iterator cell,
+  const Point<dim> &p) const
+{
+                                  // Let the start value of the
+                                  // newton iteration be the center
+                                  // of the unit cell
+  Point<dim> p_unit;
+  for (unsigned int i=0; i<dim; ++i)
+    p_unit(i)=0.5;
+
+                                  // Use the get_data function to
+                                  // create an InternalData with data
+                                  // vectors of the right size and
+                                  // transformation shape values and
+                                  // derivatives already computed at
+                                  // point p_unit.
+  const Quadrature<dim> point_quadrature(p_unit);
+  InternalData *mdata=dynamic_cast<InternalData *> (
+    get_data(update_transformation_values | update_transformation_gradients,
+            point_quadrature));
+  Assert(mdata!=0, ExcInternalError());
+
+  mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
+                                            || cell->has_boundary_lines());
+
+  MappingQ1<dim>::InternalData *p_data=0;
+  if (mdata->use_mapping_q1_on_current_cell)
+    p_data=&mdata->mapping_q1_data;
+  else
+    p_data=mdata;
+  
+                                  // perform the newton iteration.
+  transform_real_to_unit_cell_internal(cell, p, *p_data, p_unit);
+  
+  delete mdata;
+  return p_unit;
+}
+  
+  
 // explicit instantiation
 template class MappingQ<deal_II_dimension>;
index c4c31183569706d4cf5f2d65713dd117437418b2..85730cf02ed59897a775bc2b6d793e318a79619c 100644 (file)
@@ -845,15 +845,19 @@ Point<dim> MappingQ1<dim>::transform_unit_to_real_cell (
   const typename Triangulation<dim>::cell_iterator cell,
   const Point<dim> &p) const
 {
-                                  // Use the get_cell_data function
-                                  // to create an InternalData with
-                                  // data vectors already of the
-                                  // right size and with mapping
-                                  // support points set.
-  InternalData *mdata=get_cell_data(cell, update_transformation_values);
+                                  // Use the get_data function to
+                                  // create an InternalData with data
+                                  // vectors of the right size and
+                                  // transformation shape values
+                                  // already computed at point p.
+  const Quadrature<dim> point_quadrature(p);
+  InternalData *mdata=dynamic_cast<InternalData *> (
+    get_data(update_transformation_values, point_quadrature));
   Assert(mdata!=0, ExcInternalError());
 
-  compute_shapes(std::vector<Point<dim> > (1, p), *mdata);
+                                  // compute the mapping support
+                                  // points
+  compute_mapping_support_points(cell, mdata->mapping_support_points);
 
   return transform_unit_to_real_cell_internal(*mdata);
 }
@@ -876,40 +880,68 @@ Point<dim> MappingQ1<dim>::transform_unit_to_real_cell_internal (
 }
 
 
-
 template <int dim>
 Point<dim> MappingQ1<dim>::transform_real_to_unit_cell (
   const typename Triangulation<dim>::cell_iterator cell,
   const Point<dim> &p) const
 {
-                                  // Use the get_cell_data function
-                                  // to create an InternalData with
-                                  // data vectors already of the
-                                  // right size and with mapping
-                                  // support points set.
-  InternalData *mdata=get_cell_data(cell,
-                                   update_transformation_values
-                                   | update_transformation_gradients);
+                                  // Let the start value of the
+                                  // newton iteration be the center
+                                  // of the unit cell
+  Point<dim> p_unit;
+  for (unsigned int i=0; i<dim; ++i)
+    p_unit(i)=0.5;
+
+                                  // Use the get_data function to
+                                  // create an InternalData with data
+                                  // vectors of the right size and
+                                  // transformation shape values and
+                                  // derivatives already computed at
+                                  // point p_unit.
+  const Quadrature<dim> point_quadrature(p_unit);
+  InternalData *mdata=dynamic_cast<InternalData *> (
+    get_data(update_transformation_values | update_transformation_gradients,
+            point_quadrature));
   Assert(mdata!=0, ExcInternalError());
-  std::vector<Point<dim> > &points=mdata->mapping_support_points;
+
+                                  // perform the newton iteration.
+  transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
+  
+  delete mdata;
+  return p_unit;
+}
+
+
+
+template <int dim>
+void MappingQ1<dim>::transform_real_to_unit_cell_internal (
+  const typename Triangulation<dim>::cell_iterator cell,
+  const Point<dim> &p,
+  InternalData &mdata,
+  Point<dim> &p_unit) const
+{
+  const unsigned int n_shapes=mdata.shape_values.size();
+  Assert(n_shapes!=0, ExcInternalError());
+  Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
+  
+  std::vector<Point<dim> > &points=mdata.mapping_support_points;
+  compute_mapping_support_points(cell, points);
+  Assert(points.size()==n_shapes, ExcInternalError());
   
                                   // Newton iteration to solve
                                   // f(x)=p(x)-p=0
                                   // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
   
-                                  // Let the start value be the
-                                  // center of the unit cell
-                                  // (p_unit stands for x)
-  Point<dim> p_unit;
-  for (unsigned int i=0; i<dim; ++i)
-    p_unit(i)=0.5;
+                                  // The start value is set to be the
+                                  // center of the unit cell.
   
-                                  // compute shape values and
-                                  // derivatives of the mapping
-  compute_shapes(std::vector<Point<dim> > (1, p_unit), *mdata);
+                                  // The shape values and derivatives
+                                  // of the mapping at this point are
+                                  // previously computed.
+
   
                                   // f(x)
-  Point<dim> p_real(transform_unit_to_real_cell_internal(*mdata));
+  Point<dim> p_real(transform_unit_to_real_cell_internal(mdata));
   Point<dim> f = p_real-p;
 
   const double eps=1e-15*cell->diameter();
@@ -918,9 +950,9 @@ Point<dim> MappingQ1<dim>::transform_real_to_unit_cell (
     {      
                                       // f'(x)
       Tensor<2,dim> df;
-      for (unsigned int k=0; k<mdata->n_shape_functions; ++k)
+      for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
        {
-         const Tensor<1,dim> &grad_transform=mdata->derivative(0,k);
+         const Tensor<1,dim> &grad_transform=mdata.derivative(0,k);
          const Point<dim> &point=points[k];
          
          for (unsigned int i=0; i<dim; ++i)
@@ -939,40 +971,14 @@ Point<dim> MappingQ1<dim>::transform_real_to_unit_cell (
       p_unit -= d;
                                       // shape values and derivatives
                                       // at new p_unit point
-      compute_shapes(std::vector<Point<dim> > (1, p_unit), *mdata);
+      compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
       
                                       // f(x)
-      p_real=transform_unit_to_real_cell_internal(*mdata);
+      p_real=transform_unit_to_real_cell_internal(mdata);
       f = p_real-p;
     }
-
-  delete mdata;
-  return p_unit;
-}
-
-
-template <int dim>
-MappingQ1<dim>::InternalData*
-MappingQ1<dim>::get_cell_data(const typename Triangulation<dim>::cell_iterator cell,
-                             const UpdateFlags update_flags) const
-{
-  static const Point<dim> dummy_p;
-  static const Quadrature<dim> dummy_quadrature(dummy_p);
-
-  InternalData *mdata=dynamic_cast<InternalData *> (
-    get_data(update_flags, dummy_quadrature));
-  Assert(mdata!=0, ExcInternalError());
-
-                                  // compute the mapping support
-                                  // points
-  std::vector<Point<dim> > &points=mdata->mapping_support_points;
-  compute_mapping_support_points(cell, points);
-
-  return mdata;
 }
 
-  
-
 
 
 template <int dim>

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.