From 9fedb4b44d68d121d6b0d142a74b458bd802a629 Mon Sep 17 00:00:00 2001 From: heltai Date: Tue, 17 May 2011 17:26:05 +0000 Subject: [PATCH] transform_real_to_unit for codim one meshes git-svn-id: https://svn.dealii.org/trunk@23715 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 9 +++ deal.II/include/deal.II/fe/mapping.h | 9 ++- deal.II/include/deal.II/fe/mapping_q.h | 9 ++- deal.II/include/deal.II/fe/mapping_q1.h | 9 ++- deal.II/source/fe/mapping_q1.cc | 102 ++++++++++++++---------- 5 files changed, 94 insertions(+), 44 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 879fe6f756..8154b407cc 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -30,6 +30,15 @@ inconvenience this causes.
    +
  1. New: Mapping::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. +
    (Luca Heltai, 2011/05/17) + +
  2. New: The PersistentTriangulation class now works also in +the codimension one case. +
    (Luca Heltai, 2011/05/16) +
  3. Changed: Traditionally, include directories were set through the -I flag in make files in such a way that one would do @code diff --git a/deal.II/include/deal.II/fe/mapping.h b/deal.II/include/deal.II/fe/mapping.h index df8ea55093..d445dc3b8d 100644 --- a/deal.II/include/deal.II/fe/mapping.h +++ b/deal.II/include/deal.II/fe/mapping.h @@ -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 transform_real_to_unit_cell ( diff --git a/deal.II/include/deal.II/fe/mapping_q.h b/deal.II/include/deal.II/fe/mapping_q.h index 83bf78fb71..50ffaf7eed 100644 --- a/deal.II/include/deal.II/fe/mapping_q.h +++ b/deal.II/include/deal.II/fe/mapping_q.h @@ -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 * 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 transform_real_to_unit_cell ( diff --git a/deal.II/include/deal.II/fe/mapping_q1.h b/deal.II/include/deal.II/fe/mapping_q1.h index fdd07f3000..38bebfed55 100644 --- a/deal.II/include/deal.II/fe/mapping_q1.h +++ b/deal.II/include/deal.II/fe/mapping_q1.h @@ -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 * 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 transform_real_to_unit_cell ( diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index 0f541c0aa3..bacf60448e 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -1505,11 +1506,14 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it std::auto_ptr mdata(dynamic_cast ( MappingQ1::get_data(update_transformation_values - | update_transformation_gradients, - point_quadrature))); + | update_transformation_gradients, + point_quadrature))); MappingQ1::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::vertices_per_cell, ExcInternalError()); @@ -1521,35 +1525,6 @@ transform_real_to_unit_cell (const typename Triangulation::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 void MappingQ1:: @@ -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 > &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 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 &grad_transform=mdata.derivative(0,k); const Point &point=points[k]; - for (unsigned int i=0; i 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&>(f)); - - // update of p_unit - p_unit -= d; + contract (d, df_1, static_cast&>(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 > (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