--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/smartpointer.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q.h>
+
+
+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.
+ *
+ * <h3>Usage</h3>
+ *
+ * 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 <tt>dim</tt> 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<dim> fe(FE_Q<dim>(2), dim, FE_Q<dim>(1), 1);
+ * DoFHandler<dim> dof_handler(triangulation);
+ * dof_handler.distribute_dofs(fe);
+ * Vector<double> soln_vector(dof_handler.n_dofs());
+ * MappingQEulerian<dim> q2_mapping(2,soln_vector,dof_handler);
+ * @endverbatim
+ *
+ * In this example, our element consists of <tt>(dim+1)</tt> components.
+ * Only the first <tt>dim</tt> 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 <tt>EulerVectorType</tt>
+ * 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 <int dim, class EulerVectorType = Vector<double> >
+class MappingQEulerian : public MappingQ<dim>
+{
+ public:
+ /**
+ * Constructor. The first argument is
+ * the polynomical degree of the desired
+ * Qp mapping. It then takes a
+ * <tt>Vector<double> &</tt> 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
+ * <tt>DoFObjectAccessor::set_dof_values()</tt>.
+ */
+
+ MappingQEulerian (const unsigned int degree,
+ const EulerVectorType &euler_vector,
+ const DoFHandler<dim> &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<const DoFHandler<dim> > euler_dof_handler;
+
+
+ private:
+
+ /**
+ * Special quadrature rule used
+ * to define the support points
+ * in the reference configuration.
+ */
+
+ class SupportQuadrature : public Quadrature<dim>
+ {
+ 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<dim> fe_values;
+
+ /**
+ * Compute the positions of the
+ * support points in the current
+ * configuration
+ */
+ virtual void compute_mapping_support_points(
+ const typename Triangulation<dim>::cell_iterator &cell,
+ std::vector<Point<dim> > &a) const;
+
+};
+
+/*@}*/
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // __deal2__mapping_q_eulerian_h
+
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/utilities.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <lac/petsc_vector.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/fe_tools.h>
+#include <fe/mapping_q_eulerian.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+// .... MAPPING Q EULERIAN CONSTRUCTOR
+
+template <int dim, class EulerVectorType>
+MappingQEulerian<dim, EulerVectorType>::
+MappingQEulerian (const unsigned int degree,
+ const EulerVectorType &euler_vector,
+ const DoFHandler<dim> &euler_dof_handler)
+ :
+ MappingQ<dim>(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 <int dim, class EulerVectorType>
+MappingQEulerian<dim,EulerVectorType>::
+SupportQuadrature::
+SupportQuadrature (const unsigned int map_degree)
+ :
+ Quadrature<dim>(Utilities::fixed_power<dim>(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<dim> 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<unsigned int> renumber(n_q_points);
+ std::vector<unsigned int> dpo(dim+1, 1U);
+ for (unsigned int i=1; i<dpo.size(); ++i)
+ dpo[i]=dpo[i-1]*(map_degree-1);
+
+ FETools::lexicographic_to_hierarchic_numbering (
+ FiniteElementData<dim> (dpo, 1, map_degree), renumber);
+
+ // finally we assign the quadrature points in the
+ // required order.
+
+ for(unsigned int q=0; q<n_q_points; ++q)
+ this->quadrature_points[renumber[q]] = q_iterated.point(q);
+}
+
+
+
+// .... COMPUTE MAPPING SUPPORT POINTS
+
+template <int dim, class EulerVectorType>
+void
+MappingQEulerian<dim, EulerVectorType>::
+compute_mapping_support_points
+(const typename Triangulation<dim>::cell_iterator &cell,
+ std::vector<Point<dim> > &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<dim>::cell_iterator dof_cell
+ (const_cast<Triangulation<dim> *> (&(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<Vector<double> > shift_vector(n_support_pts,Vector<double>(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<n_support_pts; ++q)
+ {
+ a[q] = fe_values.quadrature_point(q);
+ for(unsigned int d=0; d<dim; ++d)
+ a[q](d) += shift_vector[q](d);
+ }
+}
+
+
+
+// explicit instantiation
+template class MappingQEulerian<deal_II_dimension, Vector<double> >;
+#ifdef DEAL_II_USE_PETSC
+template class MappingQEulerian<deal_II_dimension, PETScWrappers::Vector>;
+#endif
+
+DEAL_II_NAMESPACE_CLOSE
<li><em>Yaqi Wang:</em>
Some grid generation functions.
+
+<li><em>Joshua White:</em>
+ Arbitrary order Eulerian mappings.
</ul>
+
+
<h2>deal.II mailing list members</h2>
<p>
<ol>
+ <li> <p>New: The MappingQEulerian class provides an arbitrary order Eulerian
+ mapping where the nodes of a mesh are displaced by a previously computed
+ vector field.
+ <br>
+ (Joshua White 2008/02/05)
+ </p></li>
+
<li> <p>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
############################################################
# 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
############################################################
############################################################
############################################################
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_* \
--- /dev/null
+//---------------------------- 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 <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/numbers.h>
+#include <base/convergence_table.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_out.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <fe/mapping.h>
+#include <fe/mapping_q1.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/fe_system.h>
+#include <fe/fe_q.h>
+#include <fstream>
+#include <iostream>
+
+#include <fe/mapping_q_eulerian.h>
+
+
+// .... IMPOSED DISPLACEMENT
+
+template <int dim>
+class ImposedDisplacement : public Function<dim>
+{
+ public:
+ ImposedDisplacement() : Function<dim> (dim) { }
+ virtual void vector_value(const Point<dim> &p,
+ Vector<double> &value) const;
+};
+
+template <>
+void ImposedDisplacement<2>::vector_value(const Point<2> &p,
+ Vector<double> &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 <int dim>
+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<dim> triangulation;
+ DoFHandler<dim> dof_handler;
+ FESystem<dim> fe;
+
+ unsigned int degree;
+
+ ImposedDisplacement<dim> imposed_displacement;
+ Vector<double> displacements;
+};
+
+
+// .... CONSTRUCTOR
+
+template <int dim>
+MappingTest<dim>::MappingTest (unsigned int degree)
+ :
+ dof_handler (triangulation),
+ fe (FE_Q<dim>(degree),dim),
+ degree(degree)
+{ }
+
+
+// .... DESTRUCTOR
+
+template <int dim>
+MappingTest<dim>::~MappingTest ()
+{
+ dof_handler.clear ();
+}
+
+
+// .... COMPUTE AREA
+
+template <int dim>
+double MappingTest<dim>::compute_area ()
+{
+ QGauss<dim> quadrature_formula(degree+1);
+
+ MappingQEulerian<dim> mapping(degree,displacements,dof_handler);
+
+ FEValues<dim> 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<dim>::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<n_q_points; ++q) area += fe_values.JxW(q);
+ }
+
+ return area;
+}
+
+
+// .... RUN TEST
+
+template <int dim>
+void MappingTest<dim>::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<dim>(),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<double> (area));
+ table.add_value("error", static_cast<double> (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 <int dim>
+void MappingTest<dim>::explicitly_move_mesh ()
+{
+ std::vector<bool> moved (triangulation.n_vertices(),false);
+ unsigned int vpc = GeometryInfo<dim>::vertices_per_cell;
+
+ typename DoFHandler<dim>::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<dim> vertex_disp;
+ for (unsigned int d=0; d<dim; d++) {
+ vertex_disp[d] = displacements(cell->vertex_dof_index(v,d));
+ }
+ cell->vertex(v) += vertex_disp;
+ }}}
+}
+
+
+
+// .... GRAPHICAL OUTPUT
+
+template <int dim>
+void MappingTest<dim>::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<dim>(),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;
+}
+
--- /dev/null
+
+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::
projects / dealii-svn.git / commitdiff