From: bangerth Date: Tue, 5 Feb 2008 21:07:09 +0000 (+0000) Subject: Add Joshua White's contribution of a MappingQEulerian class. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=86004f34efc4b3fb012b845e502ad3dc9fd3ecd6;p=dealii-svn.git Add Joshua White's contribution of a MappingQEulerian class. git-svn-id: https://svn.dealii.org/trunk@15712 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/mapping_q_eulerian.h b/deal.II/deal.II/include/fe/mapping_q_eulerian.h new file mode 100644 index 0000000000..17f7e6110f --- /dev/null +++ b/deal.II/deal.II/include/fe/mapping_q_eulerian.h @@ -0,0 +1,216 @@ +//--------------------------------------------------------------------------- +// $Id: mapping_q.h 15711 2008-02-05 20:43:14Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#ifndef __deal2__mapping_q_eulerian_h +#define __deal2__mapping_q_eulerian_h + +#include +#include +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +/*!@addtogroup mapping */ +/*@{*/ + +/** + * This class is an extension of the MappingQ1Eulerian + * class to higher order Qp mappings. It is useful + * when one wants to calculate shape function information on + * a domain that is deforming as the computation proceeds. + * + *

Usage

+ * + * The constructor of this class takes three arguments: the polynomial + * degree of the desire Qp mapping, a reference to + * the vector that defines the mapping from the initial + * configuration to the current configuration, and a reference to the + * DoFHandler. The most common case is to use the solution + * vector for the problem under consideration as the shift vector. + * The key reqirement is that the number of components + * of the given vector field be equal to (or possibly greater than) the + * number of space dimensions. If there are more components than space + * dimensions (for example, if one is working with a coupled problem + * where there are additional solution variables), the + * first dim components are assumed to represent the displacement + * field, and the remaining components are ignored. If this assumption + * does not hold one may need to set up a separate DoFHandler on + * the triangulation and associate the desired shift vector to it. + * + * Typically, the DoFHandler operates on a finite element that + * is constructed as a system element (FESystem) from continuous FE_Q() + * objects. An example is shown below: + * @verbatim + * FESystem fe(FE_Q(2), dim, FE_Q(1), 1); + * DoFHandler dof_handler(triangulation); + * dof_handler.distribute_dofs(fe); + * Vector soln_vector(dof_handler.n_dofs()); + * MappingQEulerian q2_mapping(2,soln_vector,dof_handler); + * @endverbatim + * + * In this example, our element consists of (dim+1) components. + * Only the first dim components will be used, however, to define + * the Q2 mapping. The remaining components are ignored. + * + * Note that it is essential to call the distribute_dofs(...) function + * before constructing a mapping object. + * + * Also note that since the vector of shift values and the dof handler are + * only associated to this object at construction time, you have to + * make sure that whenever you use this object, the given objects + * still represent valid data. + * + * To enable the use of the MappingQ1Eulerian class also in the context + * of parallel codes using the PETSc wrapper classes, the type of + * the vector can be specified as template parameter EulerVectorType + * Not specifying this template argument in applications using the PETSc + * vector classes leads to the construction of a copy of the vector + * which is not acccessible afterwards! + * + * @author Joshua White, 2008 + */ +template > +class MappingQEulerian : public MappingQ +{ + public: + /** + * Constructor. The first argument is + * the polynomical degree of the desired + * Qp mapping. It then takes a + * Vector & to specify the + * transformation of the domain + * from the reference to + * the current configuration. + * The organization of the + * elements in the @p Vector + * must follow the concept how + * deal.II stores solutions that + * are associated to a + * triangulation. This is + * automatically the case if the + * @p Vector represents the + * solution of the previous step + * of a nonlinear problem. + * Alternatively, the @p Vector + * can be initialized by + * DoFObjectAccessor::set_dof_values(). + */ + + MappingQEulerian (const unsigned int degree, + const EulerVectorType &euler_vector, + const DoFHandler &euler_dof_handler); + + /** + * Exception + */ + + DeclException0 (ExcWrongNoOfComponents); + + /** + * Exception + */ + + DeclException0 (ExcInactiveCell); + + /** + * Exception + */ + + DeclException2 (ExcWrongVectorSize, int, int, + << "Vector has wrong size " << arg1 + << "-- expected size " << arg2); + + protected: + /** + * Reference to the vector of + * shifts. + */ + + const EulerVectorType &euler_vector; + + /** + * Pointer to the DoFHandler to + * which the mapping vector is + * associated. + */ + + const SmartPointer > euler_dof_handler; + + + private: + + /** + * Special quadrature rule used + * to define the support points + * in the reference configuration. + */ + + class SupportQuadrature : public Quadrature + { + public: + /** + * Constructor, with an argument + * defining the desired polynomial + * degree. + */ + + SupportQuadrature (const unsigned int map_degree); + + }; + + /** + * A member variable holding the + * quadrature points in the right + * order. + */ + const SupportQuadrature support_quadrature; + + /** + * FEValues object used to query the + * the given finite element field + * at the support points in the + * reference configuration. + * + * The variable is marked as + * mutable since we have to call + * FEValues::reinit from + * compute_mapping_support_points, + * a function that is 'const'. + */ + mutable FEValues fe_values; + + /** + * Compute the positions of the + * support points in the current + * configuration + */ + virtual void compute_mapping_support_points( + const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; + +}; + +/*@}*/ + + +DEAL_II_NAMESPACE_CLOSE + + +#endif // __deal2__mapping_q_eulerian_h + diff --git a/deal.II/deal.II/source/fe/mapping_q_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q_eulerian.cc new file mode 100644 index 0000000000..2177ab58ba --- /dev/null +++ b/deal.II/deal.II/source/fe/mapping_q_eulerian.cc @@ -0,0 +1,175 @@ +//--------------------------------------------------------------------------- +// $Id: mapping_q_eulerian.cc 14038 2006-10-23 02:46:34Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +// .... MAPPING Q EULERIAN CONSTRUCTOR + +template +MappingQEulerian:: +MappingQEulerian (const unsigned int degree, + const EulerVectorType &euler_vector, + const DoFHandler &euler_dof_handler) + : + MappingQ(degree, true), + euler_vector(euler_vector), + euler_dof_handler(&euler_dof_handler), + support_quadrature(degree), + fe_values(euler_dof_handler.get_fe(), + support_quadrature, + update_values | update_q_points) +{ } + + +// .... SUPPORT QUADRATURE CONSTRUCTOR + +template +MappingQEulerian:: +SupportQuadrature:: +SupportQuadrature (const unsigned int map_degree) + : + Quadrature(Utilities::fixed_power(map_degree+1)) +{ + // first we determine the support points + // on the unit cell in lexicographic order. + // for this purpose we can use an interated + // trapezoidal quadrature rule. + + QIterated q_iterated(QTrapez<1>(),map_degree); + const unsigned int n_q_points = q_iterated.n_quadrature_points; + + // we then need to define a renumbering + // vector that allows us to go from a + // lexicographic numbering scheme to a hierarchic + // one. this fragment is taking almost verbatim + // from the MappingQ class. + + std::vector renumber(n_q_points); + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i (dpo, 1, map_degree), renumber); + + // finally we assign the quadrature points in the + // required order. + + for(unsigned int q=0; qquadrature_points[renumber[q]] = q_iterated.point(q); +} + + + +// .... COMPUTE MAPPING SUPPORT POINTS + +template +void +MappingQEulerian:: +compute_mapping_support_points +(const typename Triangulation::cell_iterator &cell, + std::vector > &a) const +{ + + // first, basic assertion + // with respect to vector size, + + const unsigned int n_dofs = euler_dof_handler->n_dofs(); + const unsigned int vector_size = euler_vector.size(); + + Assert (n_dofs == vector_size,ExcWrongVectorSize(vector_size,n_dofs)); + + // we then transform our tria iterator + // into a dof iterator so we can + // access data not associated with + // triangulations + + typename DoFHandler::cell_iterator dof_cell + (const_cast *> (&(cell->get_triangulation())), + cell->level(), + cell->index(), + euler_dof_handler); + + Assert (dof_cell->active() == true, ExcInactiveCell()); + + // our quadrature rule is chosen + // so that each quadrature point + // corresponds to a support point + // in the undeformed configuration. + // we can then query the given + // displacement field at these points + // to determine the shift vector that + // maps the support points to the + // deformed configuration. + + // we assume that the given field contains + // dim displacement components, but + // that there may be other solution + // components as well (e.g. pressures). + // this class therefore assumes that the + // first dim components represent the + // actual shift vector we need, and simply + // ignores any components after that. + // this implies that the user should order + // components appropriately, or create a + // separate dof handler for the displacements. + + const unsigned int n_support_pts = support_quadrature.n_quadrature_points; + const unsigned int n_components = euler_dof_handler->get_fe().n_components(); + + Assert (n_components >= dim, ExcWrongNoOfComponents() ); + + std::vector > shift_vector(n_support_pts,Vector(n_components)); + + // fill shift vector for each support + // point using an fe_values object. + + fe_values.reinit(dof_cell); + fe_values.get_function_values(euler_vector,shift_vector); + + // and finally compute the positions of the + // support points in the deformed + // configuration. + + a.resize(n_support_pts); + for(unsigned int q=0; q >; +#ifdef DEAL_II_USE_PETSC +template class MappingQEulerian; +#endif + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/doc/authors.html b/deal.II/doc/authors.html index 8cde514ae2..1e797eedfe 100644 --- a/deal.II/doc/authors.html +++ b/deal.II/doc/authors.html @@ -134,8 +134,13 @@ alphabetical order):
  • Yaqi Wang: Some grid generation functions. + +
  • Joshua White: + Arbitrary order Eulerian mappings. + +

    deal.II mailing list members

    diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0fb045ad8d..3ebaeb842b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -233,6 +233,13 @@ constraints individually.

      +
    1. New: The MappingQEulerian class provides an arbitrary order Eulerian + mapping where the nodes of a mesh are displaced by a previously computed + vector field. +
      + (Joshua White 2008/02/05) +

    2. +
    3. New: The MappingQ class constructor now takes a parameter that determines whether a higher order mapping should also be used for interior cells. By default, the higher order mapping is only used for cells on the boundary; cells in the diff --git a/tests/fe/Makefile b/tests/fe/Makefile index 5067740a58..912cf3279a 100644 --- a/tests/fe/Makefile +++ b/tests/fe/Makefile @@ -1,6 +1,6 @@ ############################################################ # Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp -# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors ############################################################ ############################################################ @@ -22,7 +22,7 @@ default: run-tests ############################################################ tests_x=fe_data_test traits fe_tools fe_tools_* mapping \ - mapping_c1 shapes derivatives function numbering mapping_q1_eulerian \ + mapping_c1 shapes derivatives function numbering mapping_*_eulerian \ transfer internals derivatives_face \ dgq_1 \ q_* \ diff --git a/tests/fe/mapping_q_eulerian.cc b/tests/fe/mapping_q_eulerian.cc new file mode 100644 index 0000000000..cba0658aa5 --- /dev/null +++ b/tests/fe/mapping_q_eulerian.cc @@ -0,0 +1,257 @@ +//---------------------------- mapping_q_eulerian.cc --------------------------- +// mapping_q_eulerian.cc,v 1.3 2003/06/09 16:00:38 wolf Exp +// Version: +// +// Copyright (C) 2008 by the deal.II authors and Joshua White +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- mapping_q_eulerian.cc --------------------------- + +// compute some convergence results from computing pi on a mesh that +// is deformed to represent a quarter of a ring + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +// .... IMPOSED DISPLACEMENT + +template +class ImposedDisplacement : public Function +{ + public: + ImposedDisplacement() : Function (dim) { } + virtual void vector_value(const Point &p, + Vector &value) const; +}; + +template <> +void ImposedDisplacement<2>::vector_value(const Point<2> &p, + Vector &value) const +{ + double radius = 1 + (sqrt(5)-1)*p(0); + double angle = 0.5*numbers::PI*(1-p(1)); + value(0) = radius*sin(angle)-p(0); + value(1) = radius*cos(angle)-p(1); +} + + +// .... MAPPING TEST CLASS + +template +class MappingTest +{ + public: + MappingTest (unsigned int degree); + ~MappingTest (); + + void run_test(); + void graphical_output(); + + private: + double compute_area(); + void explicitly_move_mesh(); + void write_tria_to_eps(std::string id); + + Triangulation triangulation; + DoFHandler dof_handler; + FESystem fe; + + unsigned int degree; + + ImposedDisplacement imposed_displacement; + Vector displacements; +}; + + +// .... CONSTRUCTOR + +template +MappingTest::MappingTest (unsigned int degree) + : + dof_handler (triangulation), + fe (FE_Q(degree),dim), + degree(degree) +{ } + + +// .... DESTRUCTOR + +template +MappingTest::~MappingTest () +{ + dof_handler.clear (); +} + + +// .... COMPUTE AREA + +template +double MappingTest::compute_area () +{ + QGauss quadrature_formula(degree+1); + + MappingQEulerian mapping(degree,displacements,dof_handler); + + FEValues fe_values (mapping, fe, quadrature_formula, + update_JxW_values); + + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + long double area = 0.; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) { + fe_values.reinit (cell); + for(unsigned int q=0; q +void MappingTest::run_test () +{ + GridGenerator::hyper_cube (triangulation,0, 1); + + ConvergenceTable table; + + for(unsigned int ref_level = 0; + ref_level < 5; + ++ref_level, triangulation.refine_global(1)){ + + dof_handler.distribute_dofs (fe); + displacements.reinit (dof_handler.n_dofs()); + + VectorTools::interpolate(MappingQ1(),dof_handler, + imposed_displacement,displacements); + + + table.add_value("cells",triangulation.n_active_cells()); + table.add_value("dofs",dof_handler.n_dofs()); + + long double area = compute_area(); + long double error = std::fabs(numbers::PI-area)/numbers::PI; + + table.add_value("area", static_cast (area)); + table.add_value("error", static_cast (error)); + } + + table.set_precision("area", 8); + table.set_precision("error", 4); + table.set_scientific("error", true); + table.evaluate_convergence_rates("error", + ConvergenceTable::reduction_rate_log2); + table.write_text(deallog.get_file_stream()); + deallog << std::endl; + +} + + +// .... EXPLICITLY MOVE MESH + +template +void MappingTest::explicitly_move_mesh () +{ + std::vector moved (triangulation.n_vertices(),false); + unsigned int vpc = GeometryInfo::vertices_per_cell; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end(); + + for (; cell != endc; cell++) { + for (unsigned int v=0; v < vpc; v++) { + if (moved[cell->vertex_index(v)] == false){ + moved[cell->vertex_index(v)] = true; + Point vertex_disp; + for (unsigned int d=0; dvertex_dof_index(v,d)); + } + cell->vertex(v) += vertex_disp; + }}} +} + + + +// .... GRAPHICAL OUTPUT + +template +void MappingTest::graphical_output () +{ + GridGenerator::hyper_cube (triangulation,0, 1); + triangulation.refine_global(4); + + dof_handler.distribute_dofs (fe); + displacements.reinit (dof_handler.n_dofs()); + + VectorTools::interpolate(MappingQ1(),dof_handler, + imposed_displacement,displacements); + + explicitly_move_mesh(); +} + + +// .... MAIN + +int main () +{ + std::ofstream logfile ("mapping_q_eulerian/output"); + logfile.precision (2); + logfile.setf(std::ios::fixed); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + // convergence studies + + for(unsigned int degree = 1; degree <=4; ++degree) + { + deallog << ".... Q" << degree << " Mapping ...." << std::endl; + MappingTest<2> test_one(degree); + test_one.run_test(); + } + + // graphical output + + MappingTest<2> test_two(1); + test_two.graphical_output(); + + return 0; +} + diff --git a/tests/fe/mapping_q_eulerian/cmp/generic b/tests/fe/mapping_q_eulerian/cmp/generic new file mode 100644 index 0000000000..55c69da6e2 --- /dev/null +++ b/tests/fe/mapping_q_eulerian/cmp/generic @@ -0,0 +1,33 @@ + +DEAL::.... Q1 Mapping .... +cells dofs area error +1 8 2.00000000 3.6338e-01 - +4 18 2.82842712 9.9684e-02 1.87 +16 50 3.06146746 2.5505e-02 1.97 +64 162 3.12144515 6.4131e-03 1.99 +256 578 3.13654849 1.6056e-03 2.00 +DEAL:: +DEAL::.... Q2 Mapping .... +cells dofs area error +1 18 3.10456950 1.1785e-02 - +4 50 3.13914757 7.7829e-04 3.92 +16 162 3.14143772 4.9318e-05 3.98 +64 578 3.14158294 3.0930e-06 4.00 +256 2178 3.14159205 1.9348e-07 4.00 +DEAL:: +DEAL::.... Q3 Mapping .... +cells dofs area error +1 32 3.14653903 1.5745e-03 - +4 98 3.14194613 1.1251e-04 3.81 +16 338 3.14161547 7.2623e-06 3.95 +64 1250 3.14159409 4.5753e-07 3.99 +256 4802 3.14159274 2.8653e-08 4.00 +DEAL:: +DEAL::.... Q4 Mapping .... +cells dofs area error +1 50 3.14181857 7.1913e-05 - +4 162 3.14159639 1.1900e-06 5.92 +16 578 3.14159271 1.8860e-08 5.98 +64 2178 3.14159265 2.9590e-10 5.99 +256 8450 3.14159265 4.9533e-12 5.90 +DEAL::