From 646fcf1756697eba1cf6a9cd5a7c88933f8134cf Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Thu, 12 Oct 2000 13:37:56 +0000 Subject: [PATCH] New transform_real_to_unit_cell function in FiniteElement git-svn-id: https://svn.dealii.org/trunk@3433 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 9 ++++ deal.II/deal.II/include/fe/fe_system.h | 16 +++++++ deal.II/deal.II/include/fe/q1_mapping.h | 9 ++++ deal.II/deal.II/source/fe/fe_system.cc | 9 ++++ deal.II/deal.II/source/fe/q1_mapping.cc | 55 +++++++++++++++++++++++++ 5 files changed, 98 insertions(+) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index b3d0fc17c0..ebd8319e3a 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -1086,6 +1086,15 @@ class FiniteElement : public FiniteElementBase virtual Point transform_unit_to_real_cell (const DoFHandler::cell_iterator cell, const Point &p) const = 0; + /** + * 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}. + */ + virtual Point transform_real_to_unit_cell (const DoFHandler::cell_iterator cell, + const Point &p) const = 0; + /** * Return the value of the * @p{i}th shape function of the diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 9ee12608c7..74cbf06841 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -267,6 +267,22 @@ class FESystem : public FiniteElement virtual Point transform_unit_to_real_cell (const DoFHandler::cell_iterator cell, const Point &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}. As the + * transformation mapping of each + * @p{FiniteElement} of this + * @p{FESystem} should be the + * same, this function just calls + * the @p{transform} function of + * @p{base_element(0)}. + */ + virtual Point transform_real_to_unit_cell (const DoFHandler::cell_iterator cell, + const Point &p) const; + /** * Return the value of the @p{i}th shape * function of the transformation mapping diff --git a/deal.II/deal.II/include/fe/q1_mapping.h b/deal.II/deal.II/include/fe/q1_mapping.h index ee47fc6e63..62795946a2 100644 --- a/deal.II/deal.II/include/fe/q1_mapping.h +++ b/deal.II/deal.II/include/fe/q1_mapping.h @@ -56,6 +56,15 @@ class FEQ1Mapping : public FiniteElement virtual Point transform_unit_to_real_cell (const DoFHandler::cell_iterator cell, const Point &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}. + */ + virtual Point transform_real_to_unit_cell (const DoFHandler::cell_iterator cell, + const Point &p) const; + /** * Return the value of the @p{i}th shape * function at point @p{p} on the unit cell. diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 883a6d7dfe..54979a122b 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -933,6 +933,15 @@ Point FESystem::transform_unit_to_real_cell ( }; +template +Point FESystem::transform_real_to_unit_cell ( + const DoFHandler::cell_iterator cell, + const Point &p) const +{ + return base_elements[0].first->transform_real_to_unit_cell(cell, p); +}; + + template double FESystem::shape_value_transform (const unsigned int i, const Point &p) const diff --git a/deal.II/deal.II/source/fe/q1_mapping.cc b/deal.II/deal.II/source/fe/q1_mapping.cc index 3633c5c784..e0e42046e6 100644 --- a/deal.II/deal.II/source/fe/q1_mapping.cc +++ b/deal.II/deal.II/source/fe/q1_mapping.cc @@ -14,6 +14,8 @@ #include #include +#include +#include #include #include @@ -543,6 +545,59 @@ Point FEQ1Mapping::transform_unit_to_real_cell ( } +template +Point FEQ1Mapping::transform_real_to_unit_cell ( + const DoFHandler::cell_iterator cell, + const Point &p) const +{ + // Newton iteration to solve + // f(x)=p(x)-p=0 + // x_{n+1}=x_n-[f'(x)]^{-1}f(x) + + // start value in center of unit + // cell (p_unit stands for x) + Point p_unit; + for (unsigned int i=0; i p_real(transform_unit_to_real_cell(cell, p_unit)); + Point f = p_real-p; + + double eps=1e-15*cell->diameter(); + + while (f.square()>eps*eps) + { + // f'(x) + Tensor<2,dim> df; + for (unsigned int k=0; k::vertices_per_cell; ++k) + { + Tensor<1,dim> grad_transform(shape_grad_transform (k, p_unit)); + Point &vertex=cell->vertex(k); + + + for (unsigned int i=0; i d; + Tensor<2,dim> df_1; + + df_1 = invert(df); + contract (d, df_1, f); + + p_unit -= d; + // f(x) + p_real=transform_unit_to_real_cell(cell, p_unit); + f = p_real-p; + } + + return p_unit; +} + + template void FEQ1Mapping::fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, -- 2.39.5