Removed references to update_euler_vector from tets.
(Luca Heltai, 2015/04/12)
</li>
+ <li> New: A new VectorTools::get_position_vector() function has been
+ added to the library that allows one to interpolate the Geometry of
+ a (possibly curved) triangulation to vector finite element fields
+ of at least spacedim components.
+ <br>
+ (Luca Heltai, 2015/04/11)
+ <li>
+
+ <li> New: A new MappingFEField() class has been added to the library
+ that generalizes MappingQEulerian to allow arbitrary FiniteElements.
+ <br>
+ (Marco Tezzele, Luca Heltai, 2015/04/06)
+ <li>
+
<li> Changed: The cells of coarse meshes in
parallel::distributed::Triangulation used to be ordered in a Cuthill-McKee
numbering, which yields very high surface-to-volume ratios of the individual
(Wolfgang Bangerth, 2015/04/08)
</li>
- <li> New: A new MappingFEField() class has been added to the library
- that generalizes MappingQEulerian to allow arbitrary FiniteElements.
- <br>
- (Marco Tezzele, Luca Heltai, 2015/04/06)
- <li>
-
<li> Changed: In the spirit of the changes made to the distinction
between Point and Tensor objects discussed above, the first argument
to GridTools::shift() has been changed from a Point to a Tensor@<1,dim@>.
* components of the FiniteElement to use for the mapping.
*
* Typically, the DoFHandler operates on a finite element that is
- * constructed as a system element (FESystem) from continuous FE_Q()
+ * constructed as a system element (FESystem()) from continuous FE_Q()
* (for iso-parametric discretizations) or FE_Bernstein() (for
* iso-geometric discretizations) objects. An example is shown below:
*
* DoFHandler<dim,spacedim> dhq(triangulation);
* dhq.distribute_dofs(fesystem);
* Vector<double> eulerq(dhq.n_dofs());
+ * // Fills the euler vector with information from the Triangulation
+ * VectorTools::get_position_vector(dhq, eulerq);
* const ComponentMask mask(spacedim, true);
* MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
- * map.update_euler_vector_using_triangulation(eulerq);
* @endcode
*
* @author Luca Heltai, Marco Tezzele 2013, 2015
* euler_vector actually represents a valid geometry (i.e., one with
* no inverted cells, or with no zero-volume cells).
*
- * In order to facilitate the user, we provide an
- * update_euler_vector_using_triangulation() method to initialize a
- * euler_vector in a way similar to what happens in
- * MappingQEulerian, by creating a euler_vector which interpolates
- * the underlying Triangulation.
- *
* If the underlying FiniteElement is a system of FE_Q(), and
* euler_vector is initialized using
- * update_euler_vector_using_triangulation(), then this class is in
- * all respects identical to MappingQ().
+ * VectorTools::get_position_vector(), then this class is in all
+ * respects identical to MappingQ().
*
* The optional ComponentMask argument can be used to specify what
* components of the FiniteElement to use for the geometrical
*/
MappingFEField (const MappingFEField<dim,spacedim,VECTOR,DH> &mapping);
- /**
- * Helper function to fill the given euler vector with information
- * coming from the triangulation, by interpolating the geometry of
- * the Triangulation associated with the DoFHandler used at
- * construction time.
- *
- * The resulting map is guaranteed to be interpolatory at the
- * vertices of the Triangulation. Notice that this may or may not be
- * meaningful, depending on the FiniteElement you have used to
- * construct this MappingFEField.
- *
- * If the underlying FiniteElement is a system of FE_Q(), and
- * euler_vector is initialized using this function, then this class
- * is in all respects identical to MappingQ().
- */
- void update_euler_vector_using_triangulation(VECTOR &vector) const;
-
-
/**
* Transforms the point @p p on the unit cell to the point @p p_real on the
const InVector &v,
const unsigned int component);
//@}
+ /**
+ * Geometrical interpolation
+ */
+ //@{
+ /**
+ * Given a DoFHandler containing at least a spacedim vector field,
+ * this function interpolates the Triangulation at the support
+ * points of a FE_Q() finite element of the same degree as
+ * dh->get_fe().degree.
+ *
+ * Curved manifold are respected, and the resulting VECTOR will be
+ * geometrically consistent.
+ *
+ * The resulting map is guaranteed to be interpolatory at the
+ * support points of a FE_Q() finite element of the same degree as
+ * dh->get_fe().degree. Notice that this may or may not be
+ * meaningful, depending on the FiniteElement you have distribed in
+ * dh.
+ *
+ * If the underlying finite element is an FE_Q(1)^spacedim, then the
+ * resulting VECTOR is a finite element field representation of the
+ * vertices of the Triangulation.
+ *
+ * The optional ComponentMask argument can be used to specify what
+ * components of the FiniteElement to use to describe the
+ * geometry. If no mask is specified at construction time, then a
+ * default one is used, i.e., the first spacedim components of the
+ * FiniteElement are assumed to represent the geometry of the
+ * problem.
+ *
+ * @author Luca Heltai, 2015
+ */
+ template<class DH, class VECTOR>
+ void get_position_vector(const DH &dh,
+ VECTOR &vector,
+ const ComponentMask mask=ComponentMask());
+ //@}
/**
* Exception
{
return compute_mean_value(StaticMappingQ1<dim,spacedim>::mapping, dof, quadrature, v, component);
}
+
+
+ template<class DH, class VECTOR>
+ void get_position_vector(const DH &dh,
+ VECTOR &vector,
+ const ComponentMask mask)
+ {
+ AssertDimension(vector.size(), dh.n_dofs());
+
+ const unsigned int dim=DH::dimension;
+ const unsigned int spacedim=DH::space_dimension;
+ const FiniteElement<dim, spacedim> &fe = dh.get_fe();
+
+
+ // Construct default fe_mask;
+ const ComponentMask fe_mask(mask.size() ? mask :
+ ComponentMask(fe.get_nonzero_components(0).size(), true));
+
+ AssertDimension(fe_mask.size(), fe.get_nonzero_components(0).size());
+
+ std::vector<unsigned int> fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int);
+ unsigned int size = 0;
+ for (unsigned int i=0; i<fe_mask.size(); ++i)
+ {
+ if (fe_mask[i])
+ fe_to_real[i] = size++;
+ }
+ Assert(size == spacedim,
+ ExcMessage("The Component Mask you provided is invalid. It has to select exactly spacedim entries."));
+
+
+ if ( fe.has_support_points() )
+ {
+ typename DH::active_cell_iterator cell;
+ const Quadrature<dim> quad(fe.get_unit_support_points());
+
+ MappingQ<dim,spacedim> map_q(fe.degree);
+ FEValues<dim,spacedim> fe_v(map_q, fe, quad, update_quadrature_points);
+ std::vector<types::global_dof_index> dofs(fe.dofs_per_cell);
+
+ AssertDimension(fe.dofs_per_cell, fe.get_unit_support_points().size());
+ Assert(fe.is_primitive(), ExcMessage("FE is not Primitive! This won't work."));
+
+ for (cell = dh.begin_active(); cell != dh.end(); ++cell)
+ {
+ fe_v.reinit(cell);
+ cell->get_dof_indices(dofs);
+ const std::vector<Point<spacedim> > &points = fe_v.get_quadrature_points();
+ for (unsigned int q = 0; q < points.size(); ++q)
+ {
+ unsigned int comp = fe.system_to_component_index(q).first;
+ if (fe_mask[comp])
+ vector(dofs[q]) = points[q][fe_to_real[comp]];
+ }
+ }
+ }
+ else
+ {
+ // Construct a FiniteElement with FE_Q^spacedim, and call this
+ // function again.
+ //
+ // Once we have this, interpolate with the given finite element
+ // to get a Mapping which is interpolatory at the support points
+ // of FE_Q(fe.degree())
+ FESystem<dim,spacedim> feq(FE_Q<dim,spacedim>(fe.degree), spacedim);
+ DH dhq(dh.get_tria());
+ dhq.distribute_dofs(feq);
+ Vector<double> eulerq(dhq.n_dofs());
+ const ComponentMask maskq(spacedim, true);
+ get_position_vector(dhq, eulerq);
+
+ FullMatrix<double> transfer(fe.dofs_per_cell, feq.dofs_per_cell);
+ const std::vector<Point<dim> > &points = feq.get_unit_support_points();
+
+ // Here construct the interpolation matrix from FE_Q^spacedim to
+ // the FiniteElement used by euler_dof_handler.
+ //
+ // The interpolation matrix is then passed to the
+ // VectorTools::interpolate() function to generate
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ unsigned int comp_i = fe.system_to_component_index(i).first;
+ if (fe_mask[comp_i])
+ for (unsigned int j=0; j<points.size(); ++j)
+ {
+ if ( fe_to_real[comp_i] == feq.system_to_component_index(j).first)
+ transfer(i, j) = fe.shape_value(i, points[j]);
+ }
+ }
+ interpolate(dhq, dh, transfer, eulerq, vector);
+ }
+ }
}
DEAL_II_NAMESPACE_CLOSE
local_dofs[i] = (*euler_vector)(dof_indices[i]);
}
-
-template<int dim, int spacedim, class VECTOR, class DH>
-void MappingFEField<dim,spacedim,VECTOR,DH>::update_euler_vector_using_triangulation
-(VECTOR &vector) const
-{
- AssertDimension(vector.size(), euler_dof_handler->n_dofs());
- if ( fe->has_support_points() )
- {
- typename DH::active_cell_iterator cell;
- const Quadrature<dim> quad(fe->get_unit_support_points());
-
- MappingQ<dim,spacedim> map_q(fe->degree);
- FEValues<dim,spacedim> fe_v(map_q, *fe, quad, update_quadrature_points);
- std::vector<types::global_dof_index> dofs(fe->dofs_per_cell);
-
- AssertDimension(fe->dofs_per_cell, fe->get_unit_support_points().size());
- Assert(fe->is_primitive(), ExcMessage("FE is not Primitive! This won't work."));
-
- for (cell = euler_dof_handler->begin_active(); cell != euler_dof_handler->end(); ++cell)
- {
- fe_v.reinit(cell);
- cell->get_dof_indices(dofs);
- const std::vector<Point<spacedim> > &points = fe_v.get_quadrature_points();
- for (unsigned int q = 0; q < points.size(); ++q)
- {
- unsigned int comp = fe->system_to_component_index(q).first;
- if (fe_mask[comp])
- vector(dofs[q]) = points[q][fe_to_real[comp]];
- }
- }
- }
- else
- {
- // Construct a MappingFEField with an FEQ, and construct a
- // standard iso-parametric interpolation on the Triangulation.
- //
- // Once we have this, interpolate with the given finite element
- // to get a Mapping which is interpolatory at the support points
- // of FE_Q(fe->degree())
- FESystem<dim,spacedim> feq(FE_Q<dim,spacedim>(fe->degree), spacedim);
- DH dhq(euler_dof_handler->get_tria());
- dhq.distribute_dofs(feq);
- VECTOR eulerq(dhq.n_dofs());
- const ComponentMask maskq(spacedim, true);
- MappingFEField<dim,spacedim,VECTOR,DH> newfe(eulerq, dhq, maskq);
-
- newfe.update_euler_vector_using_triangulation(eulerq);
-
- FullMatrix<double> transfer(fe->dofs_per_cell, feq.dofs_per_cell);
- const std::vector<Point<dim> > &points = feq.get_unit_support_points();
-
- // Here construct the interpolation matrix from FE_Q^spacedim to
- // the FiniteElement used by euler_dof_handler.
- //
- // The interpolation matrix is then passed to the
- // VectorTools::interpolate() function to generate
- for (unsigned int i=0; i<fe->dofs_per_cell; ++i)
- {
- unsigned int comp_i = fe->system_to_component_index(i).first;
- if (fe_mask[comp_i])
- for (unsigned int j=0; j<points.size(); ++j)
- {
- if ( fe_to_real[comp_i] == feq.system_to_component_index(j).first)
- transfer(i, j) = fe->shape_value(i, points[j]);
- }
- }
- VectorTools::interpolate(dhq, *euler_dof_handler, transfer, eulerq, vector);
- }
-}
-
-
-
-
// explicit instantiations
#include "mapping_fe_field.inst"
const VEC&,
VEC&);
+ template
+ void get_position_vector
+ (const DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
+ VEC&,
+ const ComponentMask);
+
\}
#endif
}
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhb.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhb, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhb, mask);
- map.update_euler_vector_using_triangulation(eulerq);
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhb.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhb, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhb, mask);
-
- map.update_euler_vector_using_triangulation(eulerq);
typename Triangulation<dim, spacedim >::active_cell_iterator
#include <deal.II/fe/fe_bernstein.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
ComponentMask mask(spacedim+1, true);
mask.set(0, false);
+ VectorTools::get_position_vector(dhb, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhb, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhb.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhb, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhb, mask);
-
- map.update_euler_vector_using_triangulation(eulerq);
typename Triangulation<dim, spacedim >::active_cell_iterator
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhb.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhb, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhb, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhb.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhb, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhb, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhq.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhq, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/numerics/vector_tools.h>
#include <deal.II/grid/grid_out.h>
#include <fstream>
Vector<double> eulerq(dhq.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhq, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
ComponentMask mask(spacedim+1, true);
mask.set(0, false);
+ VectorTools::get_position_vector(dhq, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/numerics/vector_tools.h>
#include <deal.II/grid/grid_out.h>
#include <fstream>
Vector<double> eulerq(dhq.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhq, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhq.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhq, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
Vector<double> eulerq(dhq.n_dofs());
const ComponentMask mask(spacedim, true);
+ VectorTools::get_position_vector(dhq, eulerq, mask);
MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
- map.update_euler_vector_using_triangulation(eulerq);
-
-
typename Triangulation<dim, spacedim >::active_cell_iterator
cell = triangulation.begin_active();
projects / dealii.git / commitdiff