From: Timo Heister Date: Sat, 26 Mar 2016 05:58:29 +0000 (+0100) Subject: implement TensorProductManifold X-Git-Tag: v8.5.0-rc1~1111^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F2416%2Fhead;p=dealii.git implement TensorProductManifold --- diff --git a/include/deal.II/grid/tensor_product_manifold.h b/include/deal.II/grid/tensor_product_manifold.h new file mode 100644 index 0000000000..152c0929fe --- /dev/null +++ b/include/deal.II/grid/tensor_product_manifold.h @@ -0,0 +1,244 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__tensor_product_manifold_h +#define dealii__tensor_product_manifold_h + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +/** + * @brief Tensor product manifold of two ChartManifolds. + * + * This manifold will combine the ChartManifolds @p A and @p B given in the + * constructor to form a new ChartManifold by building the tensor product + * $A\cross B$. The first @p spacedim_A + * dimensions in the real space and the first @p chartdim_A dimensions + * of the chart will be given by manifold @p A, while the remaining + * coordinates are given by @p B. The manifold is to be used by a + * Triangulation@. + * + * An example usage would be the combination of a SphericalManifold with + * space dimension 2 and a FlatManifold with space dimension 1 to form + * a cylindrical manifold. + * + * pull_back(), push_forward(), and push_forward_gradient() are implemented + * by splitting the input argument into inputs for @p A and @p B according + * to the given dimensions and applying the corresponding operations before + * concatenating the result. + * + * @note The dimension arguments @p dim_A and @p dim_B are not used. + * + * @tparam dim Dimension of cells (needs to match first template argument of + * the Triangulation to be attached to. + * @tparam dim_A Dimension of ChartManifold A. + * @tparam spacedim_A Spacial dimension of ChartManifold A. + * @tparam chartdim_A Chart dimension of ChartManifold A. + * @tparam dim_B Dimension of ChartManifold B. + * @tparam spacedim_B Spacial dimension of ChartManifold B. + * @tparam chartdim_B Chart dimension of ChartManifold B. + * + * @author Luca Heltai, Timo Heister, 2016 + */ +template +class TensorProductManifold : + public ChartManifold +{ +public: + /** + * The chart dimension is the sum of the chart dimensions of the manifolds + * @p A and @p B. + */ + static const unsigned int chartdim = chartdim_A+chartdim_B; + /** + * The space dimension is the sum of the space dimensions of the manifolds + * @p A and @p B. + */ + static const unsigned int spacedim = spacedim_A+spacedim_B; + + /** + * Constructor. + */ + TensorProductManifold ( + const ChartManifold &manifold_A, + const ChartManifold &manifold_B); + + /** + * Pull back operation. + */ + virtual + Point + pull_back(const Point &space_point) const; + + /** + * Push forward operation. + */ + virtual + Point + push_forward(const Point &chart_point) const; + + /** + * Gradient. + */ + virtual + DerivativeForm<1,chartdim,spacedim> + push_forward_gradient(const Point &chart_point) const; + +private: + SmartPointer, + TensorProductManifold > manifold_A; + + SmartPointer, + TensorProductManifold > manifold_B; +}; + + + +/*------------------Template Implementations------------------------*/ + + + +namespace internal +{ + namespace TensorProductManifold + { + template + Tensor<1,dim1+dim2> concat(const Tensor<1,dim1> &p1, const Tensor<1,dim2> &p2) + { + Tensor<1,dim1+dim2> r; + for (unsigned int d=0; d + Point concat(const Point &p1, const Point &p2) + { + Point r; + for (unsigned int d=0; d + void split_point(const Point &source, Point &p1, Point &p2) + { + for (unsigned int d=0; d +TensorProductManifold +::TensorProductManifold( + const ChartManifold &manifold_A, + const ChartManifold &manifold_B) + : ChartManifold ( + internal::TensorProductManifold::concat( + manifold_A.get_periodicity(), + manifold_B.get_periodicity())), + manifold_A (&manifold_A), + manifold_B (&manifold_B) +{} + +template +Point::chartdim> +TensorProductManifold +::pull_back(const Point::spacedim> &space_point) const +{ + Point space_point_A; + Point space_point_B; + internal::TensorProductManifold::split_point(space_point, space_point_A, space_point_B); + + Point result_A = manifold_A->pull_back(space_point_A); + Point result_B = manifold_B->pull_back(space_point_B); + + return internal::TensorProductManifold::concat(result_A, result_B); +} + +template +Point::spacedim> +TensorProductManifold +::push_forward(const Point::chartdim> &chart_point) const +{ + Point chart_point_A; + Point chart_point_B; + internal::TensorProductManifold::split_point(chart_point, chart_point_A, chart_point_B); + + Point result_A = manifold_A->push_forward(chart_point_A); + Point result_B = manifold_B->push_forward(chart_point_B); + + return internal::TensorProductManifold::concat(result_A, result_B); +} + +template +DerivativeForm<1, + TensorProductManifold::chartdim, + TensorProductManifold::spacedim> + + TensorProductManifold + ::push_forward_gradient(const Point::chartdim> &chart_point) const +{ + Point chart_point_A; + Point chart_point_B; + internal::TensorProductManifold::split_point(chart_point, chart_point_A, chart_point_B); + + DerivativeForm<1,chartdim_A,spacedim_A> result_A + = manifold_A->push_forward_gradient(chart_point_A); + DerivativeForm<1,chartdim_B,spacedim_B> result_B + = manifold_B->push_forward_gradient(chart_point_B); + + + DerivativeForm<1,chartdim,spacedim> result; + for (unsigned int i = 0; i +#include + +#include +#include + + +void test1() +{ + const int dim=2, spacedim=2+1; + + FunctionManifold<1,2,1> F("x;x^2", "x"); + FunctionManifold<1,1,1> G("1.0+2*x", "0.5*(x-1.0)"); + + TensorProductManifold<1, 1,2,1, 1,1,1> manifold(F, G); + + // Chart points. + Point<2> cp[2]; + cp[1][0] = 1.0; + cp[1][1] = 0.5; + + // Spacedim points + std::vector > sp(2); + + // Weights + std::vector w(2); + + sp[0] = manifold.push_forward(cp[0]); + sp[1] = manifold.push_forward(cp[1]); + + for (unsigned int d=0; d<2; ++d) + if (cp[d].distance(manifold.pull_back(sp[d])) > 1e-10) + deallog << "Error! " + << cp[d] << "->" << sp[d] << "->" << manifold.pull_back(sp[d]) + << std::endl; + + unsigned int n_intermediates = 8; + + deallog << "P0: " << sp[0] + << ", P1: " << sp[1] << std::endl; + + for (unsigned int i=0; i ip = manifold.get_new_point(Quadrature(sp, w)); + Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]); + Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]); + + deallog << "P: " << ip + << ", T(P, P0): " << t1 + << ", T(P, P1): " << t2 << std::endl; + + } +} + + + +int main () +{ + initlog(); + + test1(); + + + return 0; +} + diff --git a/tests/manifold/tensor_product_manifold_01.output b/tests/manifold/tensor_product_manifold_01.output new file mode 100644 index 0000000000..83d9255c81 --- /dev/null +++ b/tests/manifold/tensor_product_manifold_01.output @@ -0,0 +1,11 @@ + +DEAL::P0: 0.00000 0.00000 1.00000, P1: 1.00000 1.00000 2.00000 +DEAL::P: 0.00000 0.00000 1.00000, T(P, P0): 0.00000 0.00000 0.00000, T(P, P1): 1.00000 0.00000 1.00000 +DEAL::P: 0.125000 0.0156250 1.12500, T(P, P0): -0.125000 -0.0312500 -0.125000, T(P, P1): 0.875000 0.218750 0.875000 +DEAL::P: 0.250000 0.0625000 1.25000, T(P, P0): -0.250000 -0.125000 -0.250000, T(P, P1): 0.750000 0.375000 0.750000 +DEAL::P: 0.375000 0.140625 1.37500, T(P, P0): -0.375000 -0.281250 -0.375000, T(P, P1): 0.625000 0.468750 0.625000 +DEAL::P: 0.500000 0.250000 1.50000, T(P, P0): -0.500000 -0.500000 -0.500000, T(P, P1): 0.500000 0.500000 0.500000 +DEAL::P: 0.625000 0.390625 1.62500, T(P, P0): -0.625000 -0.781250 -0.625000, T(P, P1): 0.375000 0.468750 0.375000 +DEAL::P: 0.750000 0.562500 1.75000, T(P, P0): -0.750000 -1.12500 -0.750000, T(P, P1): 0.250000 0.375000 0.250000 +DEAL::P: 0.875000 0.765625 1.87500, T(P, P0): -0.875000 -1.53125 -0.875000, T(P, P1): 0.125000 0.218750 0.125000 +DEAL::P: 1.00000 1.00000 2.00000, T(P, P0): -1.00000 -2.00000 -1.00000, T(P, P1): 0.00000 0.00000 0.00000 diff --git a/tests/manifold/tensor_product_manifold_02.cc b/tests/manifold/tensor_product_manifold_02.cc new file mode 100644 index 0000000000..65b79a5a18 --- /dev/null +++ b/tests/manifold/tensor_product_manifold_02.cc @@ -0,0 +1,82 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test TensorProductManifold by refining and generating normals for +// a manually constructed cylinder hull. + +#include "../tests.h" +#include +#include + + +#include +#include +#include +#include + + +void test() +{ + std::ostream &out = deallog.get_file_stream(); + + FunctionManifold<1,1> F("x","x"); + SphericalManifold<2,2> G; + + TensorProductManifold<2, 1,1,1, 2,2,2> manifold(F, G); + + // make a hull of a cylinder + Triangulation<2,3> tria; + { + Triangulation<3,3> volume_tria; + GridGenerator::cylinder(volume_tria); + std::set boundary_ids; + boundary_ids.insert(0); + GridGenerator::extract_boundary_mesh(volume_tria, tria, boundary_ids); + } + tria.set_all_manifold_ids(0); + tria.set_manifold(0, manifold); + + tria.refine_global(1); + + out << "set view equal xyz" << std::endl + << "splot '-' with lines, '-' with vectors " << std::endl; + GridOut().write_gnuplot (tria, out); + out << "e" << std::endl; + + Triangulation<2,3>::active_cell_iterator it = tria.begin_active(); + for (; it!=tria.end(); ++it) + { + Point<3> p = it->center(true); + Tensor<1,3> t1 = manifold.get_tangent_vector(p, it->vertex(0)); + Tensor<1,3> t2 = manifold.get_tangent_vector(p, it->vertex(1)); + Tensor<1,3> n = cross_product_3d(t1, t2); + n/=-n.norm(); + out << it->center() << " " << n << std::endl; + } + out << "e" << std::endl; +} + + + +int main () +{ + initlog(); + + test(); + + return 0; +} + diff --git a/tests/manifold/tensor_product_manifold_02.output b/tests/manifold/tensor_product_manifold_02.output new file mode 100644 index 0000000000..0b91e642a5 --- /dev/null +++ b/tests/manifold/tensor_product_manifold_02.output @@ -0,0 +1,261 @@ + +set view equal xyz +splot '-' with lines, '-' with vectors +-1.00000 0.707107 -0.707107 1 0 +-1.00000 -1.83697e-16 -1.00000 1 0 +-0.500000 -1.83697e-16 -1.00000 1 0 +-0.500000 0.707107 -0.707107 1 0 +-1.00000 0.707107 -0.707107 1 0 + + +-1.00000 -1.83697e-16 -1.00000 1 0 +-1.00000 -0.707107 -0.707107 1 0 +-0.500000 -0.707107 -0.707107 1 0 +-0.500000 -1.83697e-16 -1.00000 1 0 +-1.00000 -1.83697e-16 -1.00000 1 0 + + +-0.500000 0.707107 -0.707107 1 0 +-0.500000 -1.83697e-16 -1.00000 1 0 +0.00000 -1.83697e-16 -1.00000 1 0 +0.00000 0.707107 -0.707107 1 0 +-0.500000 0.707107 -0.707107 1 0 + + +-0.500000 -1.83697e-16 -1.00000 1 0 +-0.500000 -0.707107 -0.707107 1 0 +0.00000 -0.707107 -0.707107 1 0 +0.00000 -1.83697e-16 -1.00000 1 0 +-0.500000 -1.83697e-16 -1.00000 1 0 + + +-1.00000 0.707107 -0.707107 1 0 +-0.500000 0.707107 -0.707107 1 0 +-0.500000 1.00000 0.00000 1 0 +-1.00000 1.00000 0.00000 1 0 +-1.00000 0.707107 -0.707107 1 0 + + +-0.500000 0.707107 -0.707107 1 0 +0.00000 0.707107 -0.707107 1 0 +0.00000 1.00000 0.00000 1 0 +-0.500000 1.00000 0.00000 1 0 +-0.500000 0.707107 -0.707107 1 0 + + +-1.00000 1.00000 0.00000 1 0 +-0.500000 1.00000 0.00000 1 0 +-0.500000 0.707107 0.707107 1 0 +-1.00000 0.707107 0.707107 1 0 +-1.00000 1.00000 0.00000 1 0 + + +-0.500000 1.00000 0.00000 1 0 +0.00000 1.00000 0.00000 1 0 +0.00000 0.707107 0.707107 1 0 +-0.500000 0.707107 0.707107 1 0 +-0.500000 1.00000 0.00000 1 0 + + +-1.00000 -0.707107 -0.707107 1 0 +-1.00000 -1.00000 1.22465e-16 1 0 +-0.500000 -1.00000 1.22465e-16 1 0 +-0.500000 -0.707107 -0.707107 1 0 +-1.00000 -0.707107 -0.707107 1 0 + + +-1.00000 -1.00000 1.22465e-16 1 0 +-1.00000 -0.707107 0.707107 1 0 +-0.500000 -0.707107 0.707107 1 0 +-0.500000 -1.00000 1.22465e-16 1 0 +-1.00000 -1.00000 1.22465e-16 1 0 + + +-0.500000 -0.707107 -0.707107 1 0 +-0.500000 -1.00000 1.22465e-16 1 0 +0.00000 -1.00000 1.22465e-16 1 0 +0.00000 -0.707107 -0.707107 1 0 +-0.500000 -0.707107 -0.707107 1 0 + + +-0.500000 -1.00000 1.22465e-16 1 0 +-0.500000 -0.707107 0.707107 1 0 +0.00000 -0.707107 0.707107 1 0 +0.00000 -1.00000 1.22465e-16 1 0 +-0.500000 -1.00000 1.22465e-16 1 0 + + +-1.00000 0.707107 0.707107 1 0 +-0.500000 0.707107 0.707107 1 0 +-0.500000 6.12323e-17 1.00000 1 0 +-1.00000 6.12323e-17 1.00000 1 0 +-1.00000 0.707107 0.707107 1 0 + + +-0.500000 0.707107 0.707107 1 0 +0.00000 0.707107 0.707107 1 0 +0.00000 6.12323e-17 1.00000 1 0 +-0.500000 6.12323e-17 1.00000 1 0 +-0.500000 0.707107 0.707107 1 0 + + +-1.00000 6.12323e-17 1.00000 1 0 +-0.500000 6.12323e-17 1.00000 1 0 +-0.500000 -0.707107 0.707107 1 0 +-1.00000 -0.707107 0.707107 1 0 +-1.00000 6.12323e-17 1.00000 1 0 + + +-0.500000 6.12323e-17 1.00000 1 0 +0.00000 6.12323e-17 1.00000 1 0 +0.00000 -0.707107 0.707107 1 0 +-0.500000 -0.707107 0.707107 1 0 +-0.500000 6.12323e-17 1.00000 1 0 + + +0.00000 0.707107 -0.707107 1 0 +0.00000 -1.83697e-16 -1.00000 1 0 +0.500000 -1.83697e-16 -1.00000 1 0 +0.500000 0.707107 -0.707107 1 0 +0.00000 0.707107 -0.707107 1 0 + + +0.00000 -1.83697e-16 -1.00000 1 0 +0.00000 -0.707107 -0.707107 1 0 +0.500000 -0.707107 -0.707107 1 0 +0.500000 -1.83697e-16 -1.00000 1 0 +0.00000 -1.83697e-16 -1.00000 1 0 + + +0.500000 0.707107 -0.707107 1 0 +0.500000 -1.83697e-16 -1.00000 1 0 +1.00000 -1.83697e-16 -1.00000 1 0 +1.00000 0.707107 -0.707107 1 0 +0.500000 0.707107 -0.707107 1 0 + + +0.500000 -1.83697e-16 -1.00000 1 0 +0.500000 -0.707107 -0.707107 1 0 +1.00000 -0.707107 -0.707107 1 0 +1.00000 -1.83697e-16 -1.00000 1 0 +0.500000 -1.83697e-16 -1.00000 1 0 + + +0.00000 0.707107 -0.707107 1 0 +0.500000 0.707107 -0.707107 1 0 +0.500000 1.00000 0.00000 1 0 +0.00000 1.00000 0.00000 1 0 +0.00000 0.707107 -0.707107 1 0 + + +0.500000 0.707107 -0.707107 1 0 +1.00000 0.707107 -0.707107 1 0 +1.00000 1.00000 0.00000 1 0 +0.500000 1.00000 0.00000 1 0 +0.500000 0.707107 -0.707107 1 0 + + +0.00000 1.00000 0.00000 1 0 +0.500000 1.00000 0.00000 1 0 +0.500000 0.707107 0.707107 1 0 +0.00000 0.707107 0.707107 1 0 +0.00000 1.00000 0.00000 1 0 + + +0.500000 1.00000 0.00000 1 0 +1.00000 1.00000 0.00000 1 0 +1.00000 0.707107 0.707107 1 0 +0.500000 0.707107 0.707107 1 0 +0.500000 1.00000 0.00000 1 0 + + +0.00000 -0.707107 -0.707107 1 0 +0.00000 -1.00000 1.22465e-16 1 0 +0.500000 -1.00000 1.22465e-16 1 0 +0.500000 -0.707107 -0.707107 1 0 +0.00000 -0.707107 -0.707107 1 0 + + +0.00000 -1.00000 1.22465e-16 1 0 +0.00000 -0.707107 0.707107 1 0 +0.500000 -0.707107 0.707107 1 0 +0.500000 -1.00000 1.22465e-16 1 0 +0.00000 -1.00000 1.22465e-16 1 0 + + +0.500000 -0.707107 -0.707107 1 0 +0.500000 -1.00000 1.22465e-16 1 0 +1.00000 -1.00000 1.22465e-16 1 0 +1.00000 -0.707107 -0.707107 1 0 +0.500000 -0.707107 -0.707107 1 0 + + +0.500000 -1.00000 1.22465e-16 1 0 +0.500000 -0.707107 0.707107 1 0 +1.00000 -0.707107 0.707107 1 0 +1.00000 -1.00000 1.22465e-16 1 0 +0.500000 -1.00000 1.22465e-16 1 0 + + +0.00000 0.707107 0.707107 1 0 +0.500000 0.707107 0.707107 1 0 +0.500000 6.12323e-17 1.00000 1 0 +0.00000 6.12323e-17 1.00000 1 0 +0.00000 0.707107 0.707107 1 0 + + +0.500000 0.707107 0.707107 1 0 +1.00000 0.707107 0.707107 1 0 +1.00000 6.12323e-17 1.00000 1 0 +0.500000 6.12323e-17 1.00000 1 0 +0.500000 0.707107 0.707107 1 0 + + +0.00000 6.12323e-17 1.00000 1 0 +0.500000 6.12323e-17 1.00000 1 0 +0.500000 -0.707107 0.707107 1 0 +0.00000 -0.707107 0.707107 1 0 +0.00000 6.12323e-17 1.00000 1 0 + + +0.500000 6.12323e-17 1.00000 1 0 +1.00000 6.12323e-17 1.00000 1 0 +1.00000 -0.707107 0.707107 1 0 +0.500000 -0.707107 0.707107 1 0 +0.500000 6.12323e-17 1.00000 1 0 + + +e +-0.750000 0.353553 -0.853553 4.59413e-16 0.382683 -0.923880 +-0.750000 -0.353553 -0.853553 4.94753e-16 -0.382683 -0.923880 +-0.250000 0.353553 -0.853553 4.59413e-16 0.382683 -0.923880 +-0.250000 -0.353553 -0.853553 4.94753e-16 -0.382683 -0.923880 +-0.750000 0.853553 -0.353553 0.00000 0.923880 -0.382683 +-0.250000 0.853553 -0.353553 0.00000 0.923880 -0.382683 +-0.750000 0.853553 0.353553 0.00000 0.923880 0.382683 +-0.250000 0.853553 0.353553 0.00000 0.923880 0.382683 +-0.750000 -0.853553 -0.353553 -7.06790e-17 -0.923880 -0.382683 +-0.750000 -0.853553 0.353553 0.00000 -0.923880 0.382683 +-0.250000 -0.853553 -0.353553 -3.53395e-17 -0.923880 -0.382683 +-0.250000 -0.853553 0.353553 -2.47376e-16 -0.923880 0.382683 +-0.750000 0.353553 0.853553 2.47376e-16 0.382683 0.923880 +-0.250000 0.353553 0.853553 -2.47376e-16 0.382683 0.923880 +-0.750000 -0.353553 0.853553 0.00000 -0.382683 0.923880 +-0.250000 -0.353553 0.853553 0.00000 -0.382683 0.923880 +0.250000 0.353553 -0.853553 4.59413e-16 0.382683 -0.923880 +0.250000 -0.353553 -0.853553 4.94753e-16 -0.382683 -0.923880 +0.750000 0.353553 -0.853553 4.59413e-16 0.382683 -0.923880 +0.750000 -0.353553 -0.853553 4.94753e-16 -0.382683 -0.923880 +0.250000 0.853553 -0.353553 0.00000 0.923880 -0.382683 +0.750000 0.853553 -0.353553 0.00000 0.923880 -0.382683 +0.250000 0.853553 0.353553 0.00000 0.923880 0.382683 +0.750000 0.853553 0.353553 0.00000 0.923880 0.382683 +0.250000 -0.853553 -0.353553 -7.06790e-17 -0.923880 -0.382683 +0.250000 -0.853553 0.353553 0.00000 -0.923880 0.382683 +0.750000 -0.853553 -0.353553 -3.53395e-17 -0.923880 -0.382683 +0.750000 -0.853553 0.353553 -2.47376e-16 -0.923880 0.382683 +0.250000 0.353553 0.853553 2.47376e-16 0.382683 0.923880 +0.750000 0.353553 0.853553 -2.47376e-16 0.382683 0.923880 +0.250000 -0.353553 0.853553 0.00000 -0.382683 0.923880 +0.750000 -0.353553 0.853553 0.00000 -0.382683 0.923880 +e