From: Peter Munch Date: Wed, 10 Mar 2021 21:02:19 +0000 (+0100) Subject: VectorTools::get_position_vector(): add mapping as argument X-Git-Tag: v9.3.0-rc1~345^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11880%2Fhead;p=dealii.git VectorTools::get_position_vector(): add mapping as argument --- diff --git a/doc/news/changes/minor/20210310Munch b/doc/news/changes/minor/20210310Munch new file mode 100644 index 0000000000..9d08683de0 --- /dev/null +++ b/doc/news/changes/minor/20210310Munch @@ -0,0 +1,5 @@ +Improved: The function VectorTools::get_position_vector() can now also take +a Mapping object as input parameter. +
+(Peter Munch, 2021/03/10) + diff --git a/include/deal.II/numerics/vector_tools_interpolate.h b/include/deal.II/numerics/vector_tools_interpolate.h index b13b173a50..8c5b60ee58 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.h +++ b/include/deal.II/numerics/vector_tools_interpolate.h @@ -284,6 +284,21 @@ namespace VectorTools VectorType & vector, const ComponentMask & mask = ComponentMask()); + /** + * Like the above function but also taking @p mapping as argument. + * This will introduce an additional approximation between the true geometry + * specified by the manifold if the degree of the mapping is lower than the + * degree of the finite finite element in the DoFHandler @p dh, but more + * importantly it allows to fill location vectors for mappings that do not + * preserve vertex locations (like Eulerian mappings). + */ + template + void + get_position_vector(const Mapping & mapping, + const DoFHandler &dh, + VectorType & vector, + const ComponentMask & mask = ComponentMask()); + //@} } // namespace VectorTools diff --git a/include/deal.II/numerics/vector_tools_interpolate.templates.h b/include/deal.II/numerics/vector_tools_interpolate.templates.h index dfde4a80c6..c7e73c1f14 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.templates.h +++ b/include/deal.II/numerics/vector_tools_interpolate.templates.h @@ -618,11 +618,30 @@ namespace VectorTools } + template void get_position_vector(const DoFHandler &dh, VectorType & vector, const ComponentMask & mask) + { + const FiniteElement &fe = dh.get_fe(); + get_position_vector( + *fe.reference_cell().template get_default_mapping( + fe.degree), + dh, + vector, + mask); + } + + + + template + void + get_position_vector(const Mapping & map_q, + const DoFHandler &dh, + VectorType & vector, + const ComponentMask & mask) { AssertDimension(vector.size(), dh.n_dofs()); const FiniteElement &fe = dh.get_fe(); @@ -652,13 +671,7 @@ namespace VectorTools { const Quadrature quad(fe.get_unit_support_points()); - const auto map_q = - fe.reference_cell().template get_default_mapping( - fe.degree); - FEValues fe_v(*map_q, - fe, - quad, - update_quadrature_points); + FEValues fe_v(map_q, fe, quad, update_quadrature_points); std::vector dofs(fe.n_dofs_per_cell()); AssertDimension(fe.n_dofs_per_cell(), @@ -719,7 +732,7 @@ namespace VectorTools dhq.distribute_dofs(feq); Vector eulerq(dhq.n_dofs()); const ComponentMask maskq(spacedim, true); - get_position_vector(dhq, eulerq); + get_position_vector(map_q, dhq, eulerq); FullMatrix transfer(fe.n_dofs_per_cell(), feq.n_dofs_per_cell()); diff --git a/tests/vector_tools/get_position_vector_01.cc b/tests/vector_tools/get_position_vector_01.cc new file mode 100644 index 0000000000..55a02ef4e2 --- /dev/null +++ b/tests/vector_tools/get_position_vector_01.cc @@ -0,0 +1,124 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test VectorTools::get_position_vector() with Mapping as argument. + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +test() +{ + // surface mesh + const unsigned int spacedim = dim + 1; + const unsigned int fe_degree = 3; + const unsigned int n_refinements = 1; + + Triangulation tria; + GridGenerator::hyper_sphere(tria, Point(), 0.5); + tria.refine_global(n_refinements); + + FE_Q fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + FESystem fe_dim(fe, spacedim); + DoFHandler dof_handler_dim(tria); + dof_handler_dim.distribute_dofs(fe_dim); + + // set up base high-order mapping + Vector euler_vector_base(dof_handler_dim.n_dofs()); + VectorTools::get_position_vector(MappingQGeneric(4), + dof_handler_dim, + euler_vector_base); + MappingFEField mapping_base(dof_handler_dim, + euler_vector_base); + + // clear manifold + tria.reset_all_manifolds(); + + // output mesh with with MappingQGeneric(degree=4) + { + DataOutBase::VtkFlags flags; + + DataOut> data_out; + data_out.set_flags(flags); + data_out.attach_dof_handler(dof_handler); + + data_out.build_patches( + MappingQGeneric(4), + fe_degree + 1, + DataOut>::CurvedCellRegion::curved_inner_cells); + +#if false + std::ofstream output("test.0.vtk"); + data_out.write_vtk(output); +#else + data_out.write_vtk(deallog.get_file_stream()); +#endif + } + + // output mesh with with MappingFEField based on a vector constructed + // with the triangulation without manifolds and the high-order base mapping + { + Vector euler_vector(dof_handler_dim.n_dofs()); + VectorTools::get_position_vector(mapping_base, + dof_handler_dim, + euler_vector); + MappingFEField mapping(dof_handler_dim, euler_vector); + DataOutBase::VtkFlags flags; + + DataOut> data_out; + data_out.set_flags(flags); + data_out.attach_dof_handler(dof_handler); + + data_out.build_patches( + mapping, + fe_degree + 1, + DataOut>::CurvedCellRegion::curved_inner_cells); + +#if false + std::ofstream output("test.1.vtk"); + data_out.write_vtk(output); +#else + data_out.write_vtk(deallog.get_file_stream()); +#endif + } +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + test<1>(); +} diff --git a/tests/vector_tools/get_position_vector_01.output b/tests/vector_tools/get_position_vector_01.output new file mode 100644 index 0000000000..8ae4910caf --- /dev/null +++ b/tests/vector_tools/get_position_vector_01.output @@ -0,0 +1,169 @@ + +# vtk DataFile Version 3.0 +#This file was generated +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 40 double +0.353553 -0.353553 0 +0.265165 -0.390165 0 +0.176777 -0.426777 0 +0.0883883 -0.463388 0 +5.55112e-17 -0.500000 0 +5.55112e-17 -0.500000 0 +-0.0883883 -0.463388 0 +-0.176777 -0.426777 0 +-0.265165 -0.390165 0 +-0.353553 -0.353553 0 +-0.353553 -0.353553 0 +-0.390165 -0.265165 0 +-0.426777 -0.176777 0 +-0.463388 -0.0883883 0 +-0.500000 -5.55112e-17 0 +-0.500000 -5.55112e-17 0 +-0.463388 0.0883883 0 +-0.426777 0.176777 0 +-0.390165 0.265165 0 +-0.353553 0.353553 0 +0.353553 0.353553 0 +0.390165 0.265165 0 +0.426777 0.176777 0 +0.463388 0.0883883 0 +0.500000 5.55112e-17 0 +0.500000 5.55112e-17 0 +0.463388 -0.0883883 0 +0.426777 -0.176777 0 +0.390165 -0.265165 0 +0.353553 -0.353553 0 +-0.353553 0.353553 0 +-0.265165 0.390165 0 +-0.176777 0.426777 0 +-0.0883883 0.463388 0 +-5.55112e-17 0.500000 0 +-5.55112e-17 0.500000 0 +0.0883883 0.463388 0 +0.176777 0.426777 0 +0.265165 0.390165 0 +0.353553 0.353553 0 + +CELLS 32 96 +2 0 1 +2 1 2 +2 2 3 +2 3 4 +2 5 6 +2 6 7 +2 7 8 +2 8 9 +2 10 11 +2 11 12 +2 12 13 +2 13 14 +2 15 16 +2 16 17 +2 17 18 +2 18 19 +2 20 21 +2 21 22 +2 22 23 +2 23 24 +2 25 26 +2 26 27 +2 27 28 +2 28 29 +2 30 31 +2 31 32 +2 32 33 +2 33 34 +2 35 36 +2 36 37 +2 37 38 +2 38 39 + +CELL_TYPES 32 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 +POINT_DATA 40 +# vtk DataFile Version 3.0 +#This file was generated +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 40 double +0.353553 -0.353553 0 +0.277790 -0.415753 0 +0.191304 -0.461849 0 +0.0975546 -0.490409 0 +5.55112e-17 -0.500000 0 +5.55112e-17 -0.500000 0 +-0.0975546 -0.490409 0 +-0.191304 -0.461849 0 +-0.277790 -0.415753 0 +-0.353553 -0.353553 0 +-0.353553 -0.353553 0 +-0.415753 -0.277790 0 +-0.461849 -0.191304 0 +-0.490409 -0.0975546 0 +-0.500000 -5.55112e-17 0 +-0.500000 -5.55112e-17 0 +-0.490409 0.0975546 0 +-0.461849 0.191304 0 +-0.415753 0.277790 0 +-0.353553 0.353553 0 +0.353553 0.353553 0 +0.415753 0.277790 0 +0.461849 0.191304 0 +0.490409 0.0975546 0 +0.500000 5.55112e-17 0 +0.500000 5.55112e-17 0 +0.490409 -0.0975546 0 +0.461849 -0.191304 0 +0.415753 -0.277790 0 +0.353553 -0.353553 0 +-0.353553 0.353553 0 +-0.277790 0.415753 0 +-0.191304 0.461849 0 +-0.0975546 0.490409 0 +-5.55112e-17 0.500000 0 +-5.55112e-17 0.500000 0 +0.0975546 0.490409 0 +0.191304 0.461849 0 +0.277790 0.415753 0 +0.353553 0.353553 0 + +CELLS 32 96 +2 0 1 +2 1 2 +2 2 3 +2 3 4 +2 5 6 +2 6 7 +2 7 8 +2 8 9 +2 10 11 +2 11 12 +2 12 13 +2 13 14 +2 15 16 +2 16 17 +2 17 18 +2 18 19 +2 20 21 +2 21 22 +2 22 23 +2 23 24 +2 25 26 +2 26 27 +2 27 28 +2 28 29 +2 30 31 +2 31 32 +2 32 33 +2 33 34 +2 35 36 +2 36 37 +2 37 38 +2 38 39 + +CELL_TYPES 32 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 +POINT_DATA 40