//---------------------------------------------------------------------------
// $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
* 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 (
//---------------------------------------------------------------------------
// $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
* 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 (
//---------------------------------------------------------------------------
// $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
* 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 (
#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>
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());
}
-
-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>::
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.
// 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));
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);
// 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;
}
}