]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
transform_real_to_unit for codim one meshes
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 May 2011 17:26:05 +0000 (17:26 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 May 2011 17:26:05 +0000 (17:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@23715 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/mapping.h
deal.II/include/deal.II/fe/mapping_q.h
deal.II/include/deal.II/fe/mapping_q1.h
deal.II/source/fe/mapping_q1.cc

index 879fe6f7565f9f0abc24dd6bf15e4235b5d42223..8154b407cc74a6a0640461c73c089edbd26dde16 100644 (file)
@@ -30,6 +30,15 @@ inconvenience this causes.
 
 <ol>
 
+<li> New: Mapping<dim,spacedim>::transform_real_to_unit_cell  now
+works also in the codimension one case, where it performs the normal
+projection of the point on the codimension one surface.
+<br> (Luca Heltai, 2011/05/17)
+
+<li> New: The PersistentTriangulation class now works also in
+the codimension one case.
+<br> (Luca Heltai, 2011/05/16)
+
 <li> Changed: Traditionally, include directories were set through the
 <code>-I</code> flag in make files in such a way that one would do
 @code
index df8ea55093dccc18a4fa61e210635a2579c48921..d445dc3b8d4bc071df9c02295ab2bb4acb6763a6 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -245,6 +245,13 @@ class Mapping : public Subscriptor
                                      * the real cell to the point
                                      * @p p_unit on the unit cell
                                      * @p cell and returns @p p_unit.
+                                     *
+                                     * In the codimension one case,
+                                     * this function returns the
+                                     * normal projection of the real
+                                     * point @p p on the curve or
+                                     * surface identified by the @p
+                                     * cell.
                                      */
     virtual Point<dim>
     transform_real_to_unit_cell (
index 83bf78fb7128b5257f2068789310c0409d75943d..50ffaf7eeda0faf27e9b031b4534d02b469f3721 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -109,6 +109,13 @@ class MappingQ : public MappingQ1<dim,spacedim>
                                      * Uses Newton iteration and the
                                      * @p transform_unit_to_real_cell
                                      * function.
+                                     *
+                                     * In the codimension one case,
+                                     * this function returns the
+                                     * normal projection of the real
+                                     * point @p p on the curve or
+                                     * surface identified by the @p
+                                     * cell.
                                      */
     virtual Point<dim>
     transform_real_to_unit_cell (
index fdd07f3000ee8990469dab7deb0ab3c6dd5c8155..38bebfed55b73cf479b42f099fad5d5be0bb0900 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -70,6 +70,13 @@ class MappingQ1 : public Mapping<dim,spacedim>
                                      * Uses Newton iteration and the
                                      * @p transform_unit_to_real_cell
                                      * function.
+                                     *
+                                     * In the codimension one case,
+                                     * this function returns the
+                                     * normal projection of the real
+                                     * point @p p on the curve or
+                                     * surface identified by the @p
+                                     * cell.
                                      */
     virtual Point<dim>
     transform_real_to_unit_cell (
index 0f541c0aa3703ebd5d4dc5e99f2f13971d8b5524..bacf60448eca7723a6ff6e71a23aa578887c7e3c 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q1.h>
+#include <fe/mapping_q1_eulerian.h>
 
 #include <cmath>
 #include <algorithm>
@@ -1505,11 +1506,14 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
   std::auto_ptr<InternalData>
     mdata(dynamic_cast<InternalData *> (
             MappingQ1<dim,spacedim>::get_data(update_transformation_values
-                                     | update_transformation_gradients,
-                                     point_quadrature)));
+                                             | update_transformation_gradients,
+                                             point_quadrature)));
 
   MappingQ1<dim,spacedim>::compute_mapping_support_points (cell,
                                                           mdata->mapping_support_points);
+
+// compute_mapping_support_points (cell, mdata->mapping_support_points);
+  
   Assert(mdata->mapping_support_points.size() ==
          GeometryInfo<dim>::vertices_per_cell,
         ExcInternalError());
@@ -1521,35 +1525,6 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
 }
 
 
-
-template<>
-void
-MappingQ1<1,2>::
-transform_real_to_unit_cell_internal
-(const Triangulation<1,2>::cell_iterator &,
- const Point<2>                            &,
- InternalData                                     &,
- Point<1>                                       &) const
-{
-       Assert(false, ExcNotImplemented());
-}
-
-
-
-template<>
-void
-MappingQ1<2,3>::
-transform_real_to_unit_cell_internal
-(const Triangulation<2,3>::cell_iterator &,
- const Point<3>                            &,
- InternalData                                     &,
- Point<2>                                       &) const
-{
-       Assert(false, ExcNotImplemented());
-}
-
-
-
 template<int dim, int spacedim>
 void
 MappingQ1<dim,spacedim>::
@@ -1562,14 +1537,15 @@ transform_real_to_unit_cell_internal
   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<spacedim> > &points=mdata.mapping_support_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)
-
+      
                                   // The start value is set to be the
                                   // center of the unit cell.
 
@@ -1577,6 +1553,14 @@ transform_real_to_unit_cell_internal
                                   // of the mapping at this point are
                                   // previously computed.
 
+                                  // In the codimension one case, the
+                                  // f' is rectangular, and it is
+                                  // augmented with the normal to the
+                                  // cell. The tranformation from
+                                  // real to unit is intented as a
+                                  // projection to the surface (or
+                                  // curve) generated by the given
+                                  // cell.
 
                                   // f(x)
   Point<spacedim> p_real(transform_unit_to_real_cell_internal(mdata));
@@ -1587,26 +1571,55 @@ transform_real_to_unit_cell_internal
   while (f.square()>eps*eps && loop++<10)
     {
                                       // f'(x)
-      Tensor<2,dim> df;
+      Tensor<2,spacedim> df;
       for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
        {
          const Tensor<1,dim> &grad_transform=mdata.derivative(0,k);
          const Point<spacedim> &point=points[k];
 
-         for (unsigned int i=0; i<dim; ++i)
+         for (unsigned int i=0; i<spacedim; ++i)
            for (unsigned int j=0; j<dim; ++j)
              df[i][j]+=point[i]*grad_transform[j];
+         if(dim<spacedim)
+           {
+                                              // Temporarily
+                                              // transpose df, to
+                                              // access gradients as
+                                              // a single index, and
+                                              // put in the dimth row
+                                              // the normal to this
+                                              // cell.
+             df = transpose(df);
+             if ( (dim==1) && (spacedim==2) )
+               cross_product(df[1],-df[0]);
+             else if ( (dim==2) && (spacedim==3) )
+               cross_product(df[2], df[0], df[1]);
+             else
+               Assert (false, ExcNotImplemented());
+                                              // Transpose back
+             df = transpose(df);
+           }
        }
 
                                       // Solve  [f'(x)]d=f(x)
-      Tensor<1,dim> d;
-      Tensor<2,dim> df_1;
+      Tensor<1,spacedim> d;
+      Tensor<2,spacedim> df_1;
 
       df_1 = invert(df);
-      contract (d, df_1, static_cast<const Tensor<1,dim>&>(f));
-
-                                      // update of p_unit
-      p_unit -= d;
+      contract (d, df_1, static_cast<const Tensor<1,spacedim>&>(f));
+
+                                      // update of p_unit. The
+                                      // spacedimth component of
+                                      // transformed point is simply
+                                      // ignored in codimension one
+                                      // case. When this component is
+                                      // not zero, then we are
+                                      // projecting the point to the
+                                      // surface or curve identified
+                                      // by the cell.
+      for(unsigned int i=0;i<dim; ++i)
+       p_unit[i] -= d[i];
+         
                                       // shape values and derivatives
                                       // at new p_unit point
       compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
@@ -1614,6 +1627,13 @@ transform_real_to_unit_cell_internal
                                       // f(x)
       p_real = transform_unit_to_real_cell_internal(mdata);
       f = p_real-p;
+                                      // In the codimension one case,
+                                      // we only project on the
+                                      // surface, i.e., we remove the
+                                      // normal component of the
+                                      // difference.
+      if(dim<spacedim)
+       f -= (transpose(df)[dim]*f)*f;
     }
 }
 

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.