/*---------------------------- manifold.h ---------------------------*/
#include <deal.II/base/config.h>
+#include <deal.II/base/array_view.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/thread_management.h>
*/
namespace Manifolds
{
+ /**
+ * A <code>constexpr</code> helper function that returns the number of
+ * default points for the structure type pointed to by the given
+ * <code>MeshIteratorType</code>. See the documentation of
+ * Manifolds::get_default_points_and_weights() for more information.
+ */
+ template <typename MeshIteratorType>
+ inline
+ constexpr std::size_t n_default_points_per_cell()
+ {
+ // Note that in C++11 a constexpr function can only have a return
+ // statement, so we cannot alias the structure dimension
+ return GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::vertices_per_cell
+ + GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::lines_per_cell
+ + GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::quads_per_cell
+ + GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::hexes_per_cell
+ - 1; // don't count the cell itself, just the bounding objects
+ }
+
/**
* Given a general mesh iterator, construct a quadrature object that
* contains the following points:
const bool with_laplace = false) DEAL_II_DEPRECATED;
/**
- * Given a general mesh iterator, construct vectors of quadrature points and
+ * Given a general mesh iterator, construct arrays of quadrature points and
* weights that contain the following points:
* - If the iterator points to a line, then the quadrature points
* are the two vertices of the line. This results in a point vector
* <code>cell-@>face(f)</code> or <code>cell-@>line(l)</code>.
*/
template <typename MeshIteratorType>
- std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
- std::vector<double> >
+ std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
+ n_default_points_per_cell<MeshIteratorType>()>,
+ std::array<double, n_default_points_per_cell<MeshIteratorType>()> >
get_default_points_and_weights(const MeshIteratorType &iterator,
const bool with_laplace = false);
}
-
/**
* Manifolds are used to describe the geometry of boundaries of domains as
* well as the geometry of the interior. Manifold objects are therefore
*/
virtual
Point<spacedim>
- get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const;
+ get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const;
+
/**
- * Compute a new set of points that interpolate between the given points
- * @p surrounding_points. @p weights is a table with as many columns as
- * @p surrounding_points.size(). The number of rows in @p weights determines
- * how many new points will be computed and appended to the last input
- * argument @p new_points. After exit of this function, the size of
- * @p new_points equals the size at entry plus the number of rows in
- * @p weights.
+ * Compute a new set of points that interpolate between the given points @p
+ * surrounding_points. @p weights is a table with as many columns as @p
+ * surrounding_points.size(). The number of rows in @p weights must match
+ * the length of @p new_points.
*
* In its default implementation, this function simply calls get_new_point()
- * on each row of @p weights and appends those points to the output vector
+ * on each row of @p weights and writes those points into the output array
* @p new_points. However, this function is more efficient if multiple new
* points need to be generated like in MappingQGeneric and the manifold does
* expensive transformations between a chart space and the physical space,
* by implementing only the get_new_point() function.
*
* The implementation does not allow for @p surrounding_points and
- * @p new_points to point to the same vector, so make sure to pass different
+ * @p new_points to point to the same array, so make sure to pass different
* objects into the function.
*/
virtual
void
- add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const;
+ add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim>> new_points) const;
/**
* Given a point which lies close to the given manifold, it modifies it and
* the default behavior should work out of the box.
*/
virtual
- Point<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
- const Point<spacedim> &candidate) const;
+ Point<spacedim> project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const;
/**
* Backward compatibility interface. Return the point which shall become
*/
virtual
Point<spacedim>
- get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const;
+ get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const override;
+
/**
- * Compute a new set of points that interpolate between the given points
- * @p surrounding_points. @p weights is a table with as many columns as
- * @p surrounding_points.size(). The number of rows in @p weights determines
- * how many new points will be computed and appended to the last input
- * argument @p new_points. After exit of this function, the size of
- * @p new_points equals the size at entry plus the number of rows in
- * @p weights.
+ * Compute a new set of points that interpolate between the given points @p
+ * surrounding_points. @p weights is a table with as many columns as @p
+ * surrounding_points.size(). The number of rows in @p weights must match
+ * the length of @p new_points.
*
* For this particular implementation, the interpolation of the
* @p surrounding_points according to the @p weights is simply performed in
*/
virtual
void
- add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const;
+ add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim>> new_points) const;
/**
* Project to FlatManifold. This is the identity function for flat,
*/
virtual
Point<spacedim>
- project_to_manifold (const std::vector<Point<spacedim> > &points,
- const Point<spacedim> &candidate) const;
+ project_to_manifold (const ArrayView<const Point<spacedim>> &points,
+ const Point<spacedim> &candidate) const;
/**
* Return a vector that, at $\mathbf x_1$, is tangential to
*/
virtual
Point<spacedim>
- get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const;
+ get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const override;
/**
- * Compute a new set of points that interpolate between the given points
- * @p surrounding_points. @p weights is a table with as many columns as
- * @p surrounding_points.size(). The number of rows in @p weights determines
- * how many new points will be computed and appended to the last input
- * argument @p new_points. After exit of this function, the size of
- * @p new_points equals the size at entry plus the number of rows in
- * @p weights.
+ * Compute a new set of points that interpolate between the given points @p
+ * surrounding_points. @p weights is a table with as many columns as @p
+ * surrounding_points.size(). The number of rows in @p weights must match
+ * the length of @p new_points.
*
* The implementation of this function first transforms the
* @p surrounding_points to the chart space by calling pull_back(). Then, new
*/
virtual
void
- add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const;
-
+ add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim>> new_points) const;
/**
* Pull back the given point in spacedim to the Euclidean chartdim
* dimensional space.
};
-
-
/* -------------- declaration of explicit specializations ------------- */
#ifndef DOXYGEN
get_default_quadrature(const MeshIteratorType &iterator,
const bool with_laplace)
{
- const std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
- std::vector<double> > points_and_weights = get_default_points_and_weights(iterator,
- with_laplace);
- return Quadrature<MeshIteratorType::AccessorType::space_dimension>(points_and_weights.first,
- points_and_weights.second);
+ const auto points_and_weights = get_default_points_and_weights(iterator, with_laplace);
+ static const int spacedim = MeshIteratorType::AccessorType::space_dimension;
+ return Quadrature<spacedim>
+ (std::vector<Point<spacedim>>(points_and_weights.first.begin(),
+ points_and_weights.first.end()),
+ std::vector<double>(points_and_weights.second.begin(),
+ points_and_weights.second.end()));
}
+
+
template <typename MeshIteratorType>
- std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
- std::vector<double> >
+ std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
+ n_default_points_per_cell<MeshIteratorType>()>,
+ std::array<double, n_default_points_per_cell<MeshIteratorType>()> >
get_default_points_and_weights(const MeshIteratorType &iterator,
const bool with_laplace)
{
- const int spacedim = MeshIteratorType::AccessorType::space_dimension;
const int dim = MeshIteratorType::AccessorType::structure_dimension;
+ const int spacedim = MeshIteratorType::AccessorType::space_dimension;
+ constexpr std::size_t points_per_cell = n_default_points_per_cell<MeshIteratorType>();
- std::pair<std::vector<Point<spacedim> >,
- std::vector<double> > points_weights;
+ std::pair<std::array<Point<spacedim>, points_per_cell>, std::array<double, points_per_cell> >
+ points_weights;
// note that the exact weights are chosen such as to minimize the
switch (dim)
{
case 1:
- points_weights.first.resize(2);
- points_weights.second.resize(2);
+ Assert(points_weights.first.size() == 2, ExcInternalError());
+ Assert(points_weights.second.size() == 2, ExcInternalError());
points_weights.first[0] = iterator->vertex(0);
points_weights.second[0] = .5;
points_weights.first[1] = iterator->vertex(1);
points_weights.second[1] = .5;
break;
case 2:
- points_weights.first.resize(8);
- points_weights.second.resize(8);
+ Assert(points_weights.first.size() == 8, ExcInternalError());
+ Assert(points_weights.second.size() == 8, ExcInternalError());
for (unsigned int i=0; i<4; ++i)
{
GeometryInfo<dim>::vertices_per_cell+
GeometryInfo<dim>::lines_per_cell+
GeometryInfo<dim>::faces_per_cell;
- points_weights.first.resize(np);
- points_weights.second.resize(np);
- std::vector<Point<3> > *sp3 = reinterpret_cast<std::vector<Point<3> > *>(&points_weights.first);
+ Assert(points_weights.first.size() == np, ExcInternalError());
+ Assert(points_weights.second.size() == np, ExcInternalError());
+ auto *sp3 = reinterpret_cast<std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()> *>
+ (&points_weights.first);
unsigned int j=0;
#include <deal.II/base/function.h>
#include <deal.II/base/function_parser.h>
+#include <boost/container/small_vector.hpp>
+
DEAL_II_NAMESPACE_OPEN
/**
*/
virtual
Point<spacedim>
- get_new_point (const std::vector<Point<spacedim> > &vertices,
- const std::vector<double> &weights) const;
+ get_new_point (const ArrayView<const Point<spacedim>> &vertices,
+ const ArrayView<const double> &weights) const override;
/**
* The center of the spherical coordinate system.
* the base class for a detailed description of what this function does.
*/
virtual Point<spacedim>
- get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const override;
+ get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const override;
protected:
/**
- * A vector orthogonal to direcetion.
+ * A vector orthogonal to the normal direction.
*/
const Tensor<1,spacedim> normal_direction;
* the case for PolarManifold but not for Spherical manifold, so be careful
* when using the latter. In case the quality of the manifold is not good
* enough, upon mesh refinement it may happen that the transformation to a
- * chart inside the get_new_point() or add_new_points() methods produces points
- * that are outside the unit cell. Then this class throws an exception of
- * type Manifold@<dim,spacedim@>::ExcTransformationFailed. In that case,
- * the mesh should be refined before attaching this class, as done in the
- * following example:
+ * chart inside the get_new_point() or add_new_points() methods produces
+ * points that are outside the unit cell. Then this class throws an exception
+ * of type Mapping::ExcTransformationFailed. In that case, the mesh should be
+ * refined before attaching this class, as done in the following example:
*
* @code
* SphericalManifold<dim> spherical_manifold;
*/
virtual
Point<spacedim>
- get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const;
+ get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const override;
/**
- * Compute a new set of points that interpolate between the given points
- * @p surrounding_points. @p weights is a table with as many columns as
- * @p surrounding_points.size(). The number of rows in @p weights determines
- * how many new points will be computed and appended to the last input
- * argument @p new_points. After exit of this function, the size of
- * @p new_points equals the size at entry plus the number of rows in
- * @p weights.
+ * Compute a new set of points that interpolate between the given points @p
+ * surrounding_points. @p weights is a table with as many columns as @p
+ * surrounding_points.size(). The number of columns in @p weights must match
+ * the length of @p new_points.
*
* The implementation in this class overrides the method in the base class
* and computes the new point by a transfinite interpolation. The first step
*/
virtual
void
- add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const;
+ add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim>> new_points) const;
private:
/**
* 10 entries of a the indices <tt>cell->index()</tt>.
*/
std::array<unsigned int, 10>
- get_possible_cells_around_points(const std::vector<Point<spacedim> > &surrounding_points) const;
+ get_possible_cells_around_points(const ArrayView<const Point<spacedim>> &surrounding_points) const;
/**
- * Finalizes the identification of the correct chart and returns the location
- * of the surrounding points on the chart. This method internally calls
- * @p get_possible_cells_around_points().
+ * Finalizes the identification of the correct chart and populates @p
+ * chart_points with the pullbacks of the surrounding points. This method
+ * internally calls @p get_possible_cells_around_points().
+ *
+ * Returns an iterator to the cell on which the chart is defined.
*/
- std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
- std::vector<Point<dim> > >
- compute_chart_points(const std::vector<Point<spacedim> > &surrounding_points) const;
+ typename Triangulation<dim,spacedim>::cell_iterator
+ compute_chart_points(const ArrayView<const Point<spacedim>> &surrounding_points,
+ ArrayView<Point<dim>> chart_points) const;
/**
* Pull back operation into the unit coordinates on the given coarse cell.
* algorithms.
*/
virtual Point<spacedim>
- project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
- const Point<spacedim> &candidate) const;
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const;
private:
* projection algorithms.
*/
virtual Point<spacedim>
- project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
- const Point<spacedim> &candidate) const;
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const;
private:
/**
* exception is thrown.
*/
virtual Point<spacedim>
- project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
- const Point<spacedim> &candidate) const;
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const;
private:
/**
transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const Point<dim> &p) const
{
- std::vector<Point<spacedim> > vertices;
- std::vector<double> weights;
+ std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
+ std::array<double, GeometryInfo<dim>::vertices_per_cell> weights;
+
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
{
- vertices.push_back(cell->vertex(v));
- weights.push_back(GeometryInfo<dim>::d_linear_shape_function(p,v));
+ vertices[v] = cell->vertex(v);
+ weights[v] = GeometryInfo<dim>::d_linear_shape_function(p, v);
}
- return cell->get_manifold().get_new_point(vertices, weights);
+ return cell->get_manifold().get_new_point(make_array_view(vertices.begin(),
+ vertices.end()),
+ make_array_view(weights.begin(),
+ weights.end()));
}
{
for (unsigned int point=0; point<quadrature_points.size(); ++point)
{
- quadrature_points[point] = data.manifold->
- get_new_point(data.vertices,
- data.cell_manifold_quadrature_weights[point+data_set]);
+ quadrature_points[point] = data.manifold->get_new_point
+ (make_array_view(data.vertices),
+ make_array_view(data.cell_manifold_quadrature_weights[point+data_set]));
}
}
}
const Point<dim> &p = data.quad.point(point+data_set);
// And get its image on the manifold:
- const Point<spacedim> P = data.manifold->
- get_new_point(data.vertices,
- data.cell_manifold_quadrature_weights[point+data_set]);
+ const Point<spacedim> P = data.manifold->get_new_point
+ (make_array_view(data.vertices),
+ make_array_view(data.cell_manifold_quadrature_weights[point+data_set]));
// To compute the Jacobian, we choose dim points aligned
// with the dim reference axes, which are still in the
data.vertex_weights[j] = GeometryInfo<dim>::d_linear_shape_function(np, j);
const Point<spacedim> NP=
- data.manifold->get_new_point(data.vertices,
- data.vertex_weights);
+ data.manifold->get_new_point(make_array_view(data.vertices),
+ make_array_view(data.vertex_weights));
const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP);
}
else
{
- tmp_points.resize(2);
- tmp_points[0] = cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0));
- tmp_points[1] = cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 1));
- manifold.add_new_points(tmp_points,
- support_point_weights_perimeter_to_interior[0], a);
+ const std::array<Point<spacedim>, 2> vertices
+ {
+ {
+ cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
+ cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 1))
+ }
+ };
+
+ const std::size_t n_rows = support_point_weights_perimeter_to_interior[0].size(0);
+ a.resize(a.size() + n_rows);
+ auto a_view = make_array_view(a.end() - n_rows, a.end());
+ manifold.add_new_points(make_array_view(vertices.begin(),
+ vertices.end()),
+ support_point_weights_perimeter_to_interior[0],
+ a_view);
}
}
}
// extract the points surrounding a quad from the points
// already computed. First get the 4 vertices and then the points on
// the four lines
- tmp_points.resize(4 + 4*(polynomial_degree-1));
+ boost::container::small_vector<Point<3>, 200>
+ tmp_points(GeometryInfo<2>::vertices_per_cell
+ + GeometryInfo<2>::lines_per_cell*(polynomial_degree-1));
for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no,v)];
if (polynomial_degree > 1)
a[GeometryInfo<3>::vertices_per_cell +
(polynomial_degree-1)*
GeometryInfo<3>::face_to_cell_lines(face_no,line) + i];
- face->get_manifold().add_new_points (tmp_points,
+
+ const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0);
+ a.resize(a.size() + n_rows);
+ auto a_view = make_array_view(a.end() - n_rows, a.end());
+ face->get_manifold().add_new_points (make_array_view(tmp_points.begin(),
+ tmp_points.end()),
support_point_weights_perimeter_to_interior[1],
- a);
+ a_view);
}
}
}
}
else
{
- std::vector<Point<3> > vertices;
+ std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
- vertices.push_back(cell->vertex(i));
+ vertices[i] = cell->vertex(i);
+
Table<2,double> weights(Utilities::fixed_power<2>(polynomial_degree-1),
GeometryInfo<2>::vertices_per_cell);
for (unsigned int q=0, q2=0; q2<polynomial_degree-1; ++q2)
weights(q,i) = GeometryInfo<2>::d_linear_shape_function(point, i);
}
// TODO: use all surrounding points once Boundary path is removed
- cell->get_manifold().add_new_points(vertices, weights, a);
+ const std::size_t n_rows = weights.size(0);
+ a.resize(a.size() + n_rows);
+ auto a_view = make_array_view(a.end() - n_rows, a.end());
+ cell->get_manifold().add_new_points(make_array_view(vertices.begin(),
+ vertices.end()),
+ weights,
+ a_view);
}
}
if (all_manifold_ids_are_equal)
{
- std::vector<Point<spacedim> > vertices(a);
- cell->get_manifold().add_new_points(vertices, support_point_weights_cell, a);
+ const std::size_t n_rows = support_point_weights_cell.size(0);
+ a.resize(a.size() + n_rows);
+ auto a_view = make_array_view(a.end() - n_rows, a.end());
+ cell->get_manifold().add_new_points(make_array_view(a.begin(),
+ a.end() - n_rows),
+ support_point_weights_cell,
+ a_view);
}
else
switch (dim)
add_quad_support_points(cell, a);
else
{
- std::vector<Point<spacedim> > tmp_points(a);
- cell->get_manifold().add_new_points(tmp_points,
+ const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0);
+ a.resize(a.size() + n_rows);
+ auto a_view = make_array_view(a.end() - n_rows, a.end());
+ cell->get_manifold().add_new_points(make_array_view(a.begin(),
+ a.end() - n_rows),
support_point_weights_perimeter_to_interior[1],
- a);
+ a_view);
}
break;
// then compute the interior points
{
- std::vector<Point<spacedim> > tmp_points(a);
- cell->get_manifold().add_new_points(tmp_points,
+ const std::size_t n_rows = support_point_weights_perimeter_to_interior[2].size(0);
+ a.resize(a.size() + n_rows);
+ auto a_view = make_array_view(a.end() - n_rows, a.end());
+ cell->get_manifold().add_new_points(make_array_view(a.begin(),
+ a.end() - n_rows),
support_point_weights_perimeter_to_interior[2],
- a);
+ a_view);
}
break;
using namespace Manifolds;
-// This structure is used as comparison function for std::sort when sorting
-// points according to their weight.
-struct CompareWeights
-{
-public:
- CompareWeights(const std::vector<double> &weights)
- :
- compare_weights(weights)
- {}
-
- bool operator() (unsigned int a, unsigned int b) const
- {
- return compare_weights[a] < compare_weights[b];
- }
-
-private:
- const std::vector<double> &compare_weights;
-};
-
/* -------------------------- Manifold --------------------- */
-
template <int dim, int spacedim>
Manifold<dim, spacedim>::~Manifold ()
{}
template <int dim, int spacedim>
Point<spacedim>
Manifold<dim, spacedim>::
-project_to_manifold (const std::vector<Point<spacedim> > &,
+project_to_manifold (const ArrayView<const Point<spacedim>> &,
const Point<spacedim> &) const
{
Assert (false, ExcPureFunctionCalled());
const Point<spacedim> &p2,
const double w) const
{
- std::vector<Point<spacedim> > vertices;
- vertices.push_back(p1);
- vertices.push_back(p2);
- return project_to_manifold(vertices, w * p2 + (1-w)*p1 );
+ const std::array<Point<spacedim>, 2> vertices {{p1, p2}};
+ return project_to_manifold(make_array_view(vertices.begin(), vertices.end()),
+ w * p2 + (1-w)*p1);
}
template <int dim, int spacedim>
Point<spacedim>
Manifold<dim, spacedim>::
-get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const
{
const double tol = 1e-10;
const unsigned int n_points = surrounding_points.size();
// associative (as for the SphericalManifold).
boost::container::small_vector<unsigned int, 100> permutation(n_points);
std::iota(permutation.begin(), permutation.end(), 0u);
- std::sort(permutation.begin(), permutation.end(), CompareWeights(weights));
+ std::sort(permutation.begin(), permutation.end(),
+ [&weights](const std::size_t a, const std::size_t b)
+ {
+ return weights[a] < weights[b];
+ });
// Now loop over points in the order of their associated weight
Point<spacedim> p = surrounding_points[permutation[0]];
template <int dim, int spacedim>
void
Manifold<dim, spacedim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim>> new_points) const
{
AssertDimension(surrounding_points.size(), weights.size(1));
- Assert(&surrounding_points != &new_points,
- ExcMessage("surrounding_points and new_points cannot be the same "
- "array"));
- const unsigned int n_points = surrounding_points.size();
- std::vector<double> local_weights(n_points);
+ const std::size_t n_points = surrounding_points.size();
+ // TODO find a better dimension-dependent estimate for the size of this
+ // vector
+ boost::container::small_vector<double, 20> local_weights(n_points);
for (unsigned int row=0; row<weights.size(0); ++row)
{
for (unsigned int i=0; i<n_points; ++i)
local_weights[i] = weights(row,i);
- new_points.push_back(get_new_point(surrounding_points, local_weights));
+ new_points[row] = get_new_point(make_array_view(surrounding_points.begin(),
+ surrounding_points.end()),
+ make_array_view(local_weights.begin(),
+ local_weights.end()));
}
}
Manifold<dim, spacedim>::
get_new_point_on_line (const typename Triangulation<dim, spacedim>::line_iterator &line) const
{
- const std::pair<std::vector<Point<spacedim> >, std::vector<double> > points_weights(get_default_points_and_weights(line));
- return get_new_point (points_weights.first,points_weights.second);
+ const auto points_weights = get_default_points_and_weights(line);
+ return get_new_point (make_array_view(points_weights.first.begin(),
+ points_weights.first.end()),
+ make_array_view(points_weights.second.begin(),
+ points_weights.second.end()));
}
Manifold<dim, spacedim>::
get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
{
- const std::pair<std::vector<Point<spacedim> >, std::vector<double> > points_weights(get_default_points_and_weights(quad));
- return get_new_point (points_weights.first,points_weights.second);
+ const auto points_weights = get_default_points_and_weights(quad);
+ return get_new_point (make_array_view(points_weights.first.begin(),
+ points_weights.first.end()),
+ make_array_view(points_weights.second.begin(),
+ points_weights.second.end()));
}
Manifold<3,3>::
get_new_point_on_hex (const Triangulation<3, 3>::hex_iterator &hex) const
{
- const std::pair<std::vector<Point<3> >, std::vector<double> > points_weights(get_default_points_and_weights(hex,true));
- return get_new_point (points_weights.first,points_weights.second);
+ const auto points_weights = get_default_points_and_weights(hex, true);
+ return get_new_point (make_array_view(points_weights.first.begin(),
+ points_weights.first.end()),
+ make_array_view(points_weights.second.begin(),
+ points_weights.second.end()));
}
{
const double epsilon = 1e-8;
- std::vector<Point<spacedim> > q;
- q.push_back(x1);
- q.push_back(x2);
-
- std::vector<double> w;
- w.push_back(epsilon);
- w.push_back(1.0-epsilon);
-
- const Tensor<1,spacedim> neighbor_point = get_new_point (q, w);
+ const std::array<Point<spacedim>, 2> points {{x1, x2}};
+ const std::array<double, 2> weights {{epsilon, 1.0 - epsilon}};
+ const Point<spacedim> neighbor_point = get_new_point (make_array_view(points.begin(),
+ points.end()),
+ make_array_view(weights.begin(),
+ weights.end()));
return (neighbor_point-x1)/epsilon;
}
template <int dim, int spacedim>
Point<spacedim>
FlatManifold<dim, spacedim>::
-get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const
{
Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0)-1.0) < 1e-10,
ExcMessage("The weights for the individual points should sum to 1!"));
template <int dim, int spacedim>
void
FlatManifold<dim, spacedim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim>> new_points) const
{
AssertDimension(surrounding_points.size(), weights.size(1));
if (weights.size(0) == 0)
return;
- const unsigned int n_points = surrounding_points.size();
+ const std::size_t n_points = surrounding_points.size();
Tensor<1,spacedim> minP = periodicity;
for (unsigned int d=0; d<spacedim; ++d)
// check whether periodicity shifts some of the points. Only do this if
// necessary to avoid memory allocation
const Point<spacedim> *surrounding_points_start = &surrounding_points[0];
- std::vector<Point<spacedim> > modified_points;
+
+ boost::container::small_vector<Point<spacedim>, 200> modified_points;
bool adjust_periodicity = false;
for (unsigned int d=0; d<spacedim; ++d)
if (periodicity[d] > 0)
}
if (adjust_periodicity == true)
{
- modified_points = surrounding_points;
+ modified_points.resize(surrounding_points.size());
+ std::copy(surrounding_points.begin(), surrounding_points.end(),
+ modified_points.begin());
for (unsigned int d=0; d<spacedim; ++d)
if (periodicity[d] > 0)
for (unsigned int i=0; i<n_points; ++i)
if (new_point[d] < 0)
new_point[d] += periodicity[d];
- new_points.push_back(project_to_manifold(surrounding_points, new_point));
+ // TODO should this use surrounding_points_start or surrounding_points?
+ // The older version used surrounding_points
+ new_points[row] = project_to_manifold(make_array_view(surrounding_points.begin(),
+ surrounding_points.end()),
+ new_point);
}
}
-
template <int dim, int spacedim>
Point<spacedim>
-FlatManifold<dim, spacedim>::project_to_manifold (const std::vector<Point<spacedim> > &/*vertices*/,
- const Point<spacedim> &candidate) const
+FlatManifold<dim, spacedim>::project_to_manifold
+(const ArrayView<const Point<spacedim>> &/*vertices*/,
+ const Point<spacedim> &candidate) const
{
return candidate;
}
template <int dim, int spacedim, int chartdim>
Point<spacedim>
ChartManifold<dim,spacedim,chartdim>::
-get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const
{
- std::vector<Point<chartdim> > chart_points(surrounding_points.size());
+ const std::size_t n_points = surrounding_points.size();
+
+ boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
- for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ for (unsigned int i=0; i<n_points; ++i)
chart_points[i] = pull_back(surrounding_points[i]);
- const Point<chartdim> p_chart = sub_manifold.get_new_point(chart_points,weights);
+ const Point<chartdim> p_chart = sub_manifold.get_new_point
+ (make_array_view(chart_points.begin(), chart_points.end()),
+ weights);
return push_forward(p_chart);
}
template <int dim, int spacedim, int chartdim>
void
ChartManifold<dim,spacedim,chartdim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim>> new_points) const
{
Assert(weights.size(0) > 0, ExcEmptyObject());
AssertDimension(surrounding_points.size(), weights.size(1));
- const unsigned int n_points = surrounding_points.size();
+ const std::size_t n_points = surrounding_points.size();
- std::vector<Point<chartdim> > chart_points(n_points);
- for (unsigned int i=0; i<n_points; ++i)
+ boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
+ for (std::size_t i=0; i<n_points; ++i)
chart_points[i] = pull_back(surrounding_points[i]);
- std::vector<Point<chartdim> > new_points_on_chart;
- new_points_on_chart.reserve(weights.size(0));
- sub_manifold.add_new_points(chart_points, weights, new_points_on_chart);
+ boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(weights.size(0));
+ sub_manifold.add_new_points(make_array_view(chart_points.begin(),
+ chart_points.end()),
+ weights,
+ make_array_view(new_points_on_chart.begin(),
+ new_points_on_chart.end()));
- for (unsigned int row=0; row<weights.size(0); ++row)
- new_points.push_back(push_forward(new_points_on_chart[row]));
+ for (std::size_t row=0; row<weights.size(0); ++row)
+ new_points[row] = push_forward(new_points_on_chart[row]);
}
template <int dim, int spacedim>
Point<spacedim>
SphericalManifold<dim,spacedim>::
-get_new_point (const std::vector<Point<spacedim> > &vertices,
- const std::vector<double> &weights) const
+get_new_point (const ArrayView<const Point<spacedim>> &vertices,
+ const ArrayView<const double> &weights) const
{
const unsigned int n_points = vertices.size();
template <int dim, int spacedim>
Point<spacedim>
CylindricalManifold<dim,spacedim>::
-get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const
{
// First check if the average in space lies on the axis.
Point<spacedim> middle;
FunctionManifold<dim,spacedim,chartdim>::push_forward_gradient(const Point<chartdim> &chart_point) const
{
DerivativeForm<1, chartdim, spacedim> DF;
- std::vector<Tensor<1, chartdim> > gradients(spacedim);
- push_forward_function->vector_gradient(chart_point, gradients);
for (unsigned int i=0; i<spacedim; ++i)
- for (unsigned int j=0; j<chartdim; ++j)
- DF[i][j] = gradients[i][j];
+ {
+ const auto gradient = push_forward_function->gradient(chart_point, i);
+ for (unsigned int j=0; j<chartdim; ++j)
+ DF[i][j] = gradient[j];
+ }
return DF;
}
// add the contribution from the lines around the cell (first line in
// formula)
- std::vector<double> weights(GeometryInfo<2>::vertices_per_face);
- std::vector<Point<spacedim> > points(GeometryInfo<2>::vertices_per_face);
+ std::array<double, GeometryInfo<2>::vertices_per_face> weights;
+ std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_face> points;
+ // note that the views are immutable, but the arrays are not
+ const auto weights_view = make_array_view(weights.begin(), weights.end());
+ const auto points_view = make_array_view(points.begin(), points.end());
+
for (unsigned int line=0; line<GeometryInfo<2>::lines_per_cell; ++line)
{
const double my_weight = line%2 ? chart_point[line/2] : 1-chart_point[line/2];
weights[0] = 1. - line_point;
weights[1] = line_point;
new_point += my_weight *
- cell.line(line)->get_manifold().get_new_point(points, weights);
+ cell.line(line)->get_manifold().get_new_point(points_view,
+ weights_view);
}
}
weights_lines[line] = 0;
// start with the contributions of the faces
- std::vector<double> weights;
- std::vector<Point<spacedim> > points;
+ std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
+ std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> points;
+ // note that the views are immutable, but the arrays are not
+ const auto weights_view = make_array_view(weights.begin(), weights.end());
+ const auto points_view = make_array_view(points.begin(), points.end());
+
for (unsigned int face=0; face<GeometryInfo<3>::faces_per_cell; ++face)
{
Point<2> quad_point(chart_point[(face/2+1)%3], chart_point[(face/2+2)%3]);
}
else
{
- points.resize(GeometryInfo<2>::vertices_per_cell);
- weights.resize(GeometryInfo<2>::vertices_per_cell);
for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
{
points[v] = cell.vertex(GeometryInfo<3>::face_to_cell_vertices(face,v));
weights[v] = GeometryInfo<2>::d_linear_shape_function(quad_point, v);
}
new_point += my_weight *
- cell.face(face)->get_manifold().get_new_point(points, weights);
+ cell.face(face)->get_manifold().get_new_point(points_view,
+ weights_view);
}
}
}
else
{
- points.resize(2);
- weights.resize(2);
points[0] = cell.vertex(GeometryInfo<3>::line_to_cell_vertices(line,0));
points[1] = cell.vertex(GeometryInfo<3>::line_to_cell_vertices(line,1));
weights[0] = 1. - line_point;
weights[1] = line_point;
new_point -= my_weight *
- cell.line(line)->get_manifold().get_new_point(points, weights);
+ cell.line(line)->get_manifold().get_new_point(points_view,
+ weights_view);
}
}
template <int dim, int spacedim>
std::array<unsigned int, 10>
TransfiniteInterpolationManifold<dim,spacedim>
-::get_possible_cells_around_points(const std::vector<Point<spacedim> > &points) const
+::get_possible_cells_around_points(const ArrayView<const Point<spacedim>> &points) const
{
// The methods to identify cells around points in GridTools are all written
// for the active cells, but we are here looking at some cells at the coarse
typename Triangulation<dim,spacedim>::cell_iterator
cell = triangulation->begin(level_coarse),
endc = triangulation->end(level_coarse);
- std::vector<std::pair<double, unsigned int> > distances_and_cells;
+ boost::container::small_vector<std::pair<double, unsigned int>, 200> distances_and_cells;
for ( ; cell != endc; ++cell)
{
// only consider cells where the current manifold is attached
template <int dim, int spacedim>
-std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
- std::vector<Point<dim> > >
- TransfiniteInterpolationManifold<dim, spacedim>
- ::compute_chart_points (const std::vector<Point<spacedim> > &surrounding_points) const
+typename Triangulation<dim,spacedim>::cell_iterator
+TransfiniteInterpolationManifold<dim, spacedim>
+::compute_chart_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ ArrayView<Point<dim>> chart_points) const
{
- std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
- std::vector<Point<dim> > > chart_points;
- chart_points.second.resize(surrounding_points.size());
+ Assert(surrounding_points.size() == chart_points.size(),
+ ExcMessage("The chart points array view must be as large as the "
+ "surrounding points array view."));
std::array<unsigned int,10> nearby_cells =
get_possible_cells_around_points(surrounding_points);
bool inside_unit_cell = true;
for (unsigned int i=0; i<surrounding_points.size(); ++i)
{
- chart_points.second[i] = pull_back(cell, surrounding_points[i]);
+ chart_points[i] = pull_back(cell, surrounding_points[i]);
// Tolerance 1e-6 chosen that the method also works with
// SphericalManifold
- if (GeometryInfo<dim>::is_inside_unit_cell(chart_points.second[i],
+ if (GeometryInfo<dim>::is_inside_unit_cell(chart_points[i],
1e-6) == false)
{
inside_unit_cell = false;
}
if (inside_unit_cell == true)
{
- chart_points.first = cell;
- return chart_points;
+ return cell;
}
}
// a valid inversion should have returned a point above.
AssertThrow(false,
(typename Mapping<dim,spacedim>::ExcTransformationFailed()));
- chart_points.second.clear();
- return chart_points;
+ return typename Triangulation<dim,spacedim>::cell_iterator();
}
template <int dim, int spacedim>
Point<spacedim>
TransfiniteInterpolationManifold<dim, spacedim>
-::get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
- const std::vector<double> &weights) const
+::get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const ArrayView<const double> &weights) const
{
- const std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
- std::vector<Point<dim> > > chart_points =
- compute_chart_points(surrounding_points);
+ boost::container::small_vector<Point<dim>, 100> chart_points(surrounding_points.size());
+ ArrayView<Point<dim>> chart_points_view = make_array_view(chart_points.begin(),
+ chart_points.end());
+ const auto cell = compute_chart_points(surrounding_points, chart_points_view);
- const Point<dim> p_chart = chart_manifold.get_new_point(chart_points.second,weights);
+ const Point<dim> p_chart = chart_manifold.get_new_point (chart_points_view, weights);
- return push_forward(chart_points.first, p_chart);
+ return push_forward(cell, p_chart);
}
template <int dim, int spacedim>
void
TransfiniteInterpolationManifold<dim,spacedim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
- const Table<2,double> &weights,
- std::vector<Point<spacedim> > &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Table<2,double> &weights,
+ ArrayView<Point<spacedim> > new_points) const
{
Assert(weights.size(0) > 0, ExcEmptyObject());
AssertDimension(surrounding_points.size(), weights.size(1));
- const std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
- std::vector<Point<dim> > > chart_points =
- compute_chart_points(surrounding_points);
+ boost::container::small_vector<Point<dim>, 100> chart_points(surrounding_points.size());
+ ArrayView<Point<dim>> chart_points_view = make_array_view(chart_points.begin(),
+ chart_points.end());
+ const auto cell = compute_chart_points(surrounding_points, chart_points_view);
- std::vector<Point<dim> > new_points_on_chart;
- new_points_on_chart.reserve(weights.size(0));
- chart_manifold.add_new_points(chart_points.second, weights, new_points_on_chart);
+ boost::container::small_vector<Point<dim>, 100> new_points_on_chart(weights.size(0));
+ chart_manifold.add_new_points(chart_points_view,
+ weights,
+ make_array_view(new_points_on_chart.begin(),
+ new_points_on_chart.end()));
for (unsigned int row=0; row<weights.size(0); ++row)
- new_points.push_back(push_forward(chart_points.first, new_points_on_chart[row]));
+ new_points[row] = push_forward(cell, new_points_on_chart[row]);
}
// as returned by the underlying manifold
// object.
{
- std::vector<Point<spacedim> > ps(2);
- std::vector<double> ws(2, 0.5);
- ps[0] = cell->face(boundary_face)
- ->child(0)->vertex(1);
- ps[1] = cell->face(GeometryInfo<dim>
- ::opposite_face[boundary_face])
- ->child(0)->vertex(1);
+ const std::array<Point<spacedim>, 2> ps
+ {
+ {
+ cell->face(boundary_face)->child(0)->vertex(1),
+ cell->face(GeometryInfo<dim>::opposite_face[boundary_face])
+ ->child(0)->vertex(1)
+ }
+ };
+ const std::array<double, 2> ws {{0.5, 0.5}};
triangulation.vertices[next_unused_vertex]
- = cell->get_manifold().get_new_point(ps,ws);
+ = cell->get_manifold().get_new_point(make_array_view(ps.begin(),
+ ps.end()),
+ make_array_view(ws.begin(),
+ ws.end()));
}
}
}
else
{
TriaRawIterator<TriaAccessor<structdim, dim, spacedim> > it(obj);
- const std::pair<std::vector<Point<spacedim> >,
- std::vector<double> > points_and_weights = Manifolds::get_default_points_and_weights(it, use_laplace);
- return obj.get_manifold().get_new_point(points_and_weights.first,
- points_and_weights.second);
+ const auto points_and_weights = Manifolds::get_default_points_and_weights(it, use_laplace);
+ return obj.get_manifold().get_new_point(make_array_view(points_and_weights.first.begin(),
+ points_and_weights.first.end()),
+ make_array_view(points_and_weights.second.begin(),
+ points_and_weights.second.end()));
}
}
}
TriaAccessor<structdim, dim, spacedim>::intermediate_point (const Point<structdim> &coordinates) const
{
// Surrounding points and weights.
- std::vector<Point<spacedim> > p(GeometryInfo<structdim>::vertices_per_cell);
- std::vector<double> w(GeometryInfo<structdim>::vertices_per_cell);
+ std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
+ std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
{
w[i] = GeometryInfo<structdim>::d_linear_shape_function(coordinates, i);
}
- return this->get_manifold().get_new_point(p, w);
+ return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
+ make_array_view(w.begin(), w.end()));
}
template <int dim, int spacedim>
Point<spacedim> NormalProjectionBoundary<dim,spacedim>::
- project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
- const Point<spacedim> &candidate) const
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const
{
(void)surrounding_points;
#ifdef DEBUG
template <int dim, int spacedim>
Point<spacedim> DirectionalProjectionBoundary<dim,spacedim>::
- project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
- const Point<spacedim> &candidate) const
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const
{
(void)surrounding_points;
#ifdef DEBUG
template <int dim, int spacedim>
Point<spacedim> NormalToMeshProjectionBoundary<dim,spacedim>::
- project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
- const Point<spacedim> &candidate) const
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const
{
TopoDS_Shape out_shape;
Tensor<1,3> average_normal;
for (unsigned int i=0; i<ps.size(); ++i)
{
- middle = manifold.get_new_point(ps[i],ws);
+ middle = manifold.get_new_point(make_array_view(ps[i]),
+ make_array_view(ws));
deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
}
return 0;
}
-
for (unsigned int i=0; i<ps.size(); ++i)
{
- middle = manifold.get_new_point(ps[i],ws);
+ middle = manifold.get_new_point(make_array_view(ps[i]),
+ make_array_view(ws));
deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
}
return 0;
}
-
weights[0] = (double)i/((double)n_intermediates-1);
weights[1] = 1.0-weights[0];
- deallog << manifold.get_new_point(points, weights) << std::endl;
+ deallog << manifold.get_new_point(make_array_view(points),
+ make_array_view(weights))
+ << std::endl;
}
}
return 0;
}
-
-
-
w[0] = 1.0-(double)i/((double)n_intermediates);
w[1] = 1.0 - w[0];
- Point<spacedim> ip = manifold.get_new_point(sp, w);
+ Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+ make_array_view(w));
Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
return 0;
}
-
w[0] = 1.0-(double)i/((double)n_intermediates);
w[1] = 1.0 - w[0];
- Point<spacedim> ip = manifold.get_new_point(sp, w);
+ Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+ make_array_view(w));
Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
return 0;
}
-
w[0] = 1.0-(double)i/((double)n_intermediates);
w[1] = 1.0 - w[0];
- Point<spacedim> ip = manifold.get_new_point(sp, w);
+ Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+ make_array_view(w));
Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
return 0;
}
-
w[0] = 1.0-(double)i/((double)n_intermediates);
w[1] = 1.0 - w[0];
- Point<spacedim> ip = manifold.get_new_point(sp, w);
+ Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+ make_array_view(w));
Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
w[1] = 1.0 - w[0];
Point<spacedim> ip = manifold.
- pull_back(manifold.get_new_point(sp, w));
+ pull_back(manifold.get_new_point(make_array_view(sp),
+ make_array_view(w)));
ip[0] = w[1];
for (unsigned int i=0; i<ps.size(); ++i)
{
- middle = manifold.get_new_point(ps[i],ws);
+ middle = manifold.get_new_point(make_array_view(ps[i]),
+ make_array_view(ws));
deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
}
return 0;
}
-
weights[0] = (double)i/((double)n_intermediates-1);
weights[1] = 1.0-weights[0];
- deallog << manifold.get_new_point(points, weights) << std::endl;
+ deallog << manifold.get_new_point(make_array_view(points),
+ make_array_view(weights))
+ << std::endl;
}
}
return 0;
}
-
-
-
w[0] = 1.0-(double)i/((double)n_intermediates);
w[1] = 1.0 - w[0];
- Point<spacedim> ip = manifold.get_new_point(p, w);
+ Point<spacedim> ip = manifold.get_new_point(make_array_view(p),
+ make_array_view(w));
Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, p[0]);
Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, p[1]);
return 0;
}
-
weights[0] = (double)i/((double)n_intermediates-1);
weights[1] = 1.0-weights[0];
- deallog << manifold.get_new_point(points, weights) << std::endl;
+ deallog << manifold.get_new_point(make_array_view(points),
+ make_array_view(weights))
+ << std::endl;
}
}
return 0;
}
-
-
-
{
public:
Point<spacedim>
- project_to_manifold (const std::vector<Point<spacedim> > &vertices,
- const Point<spacedim> &candidate) const
+ project_to_manifold (const ArrayView<const Point<spacedim> > &vertices,
+ const Point<spacedim> &candidate) const override
{
// Shift the y coordinate to 4*x*(1-x)
Point<spacedim> p = candidate;
return 0;
}
-
weights[1] = 1.0/3.0;
weights[2] = 1.0/3.0;
- Point<3> Q = manifold.get_new_point(points1, weights);
- Point<3> S = manifold.get_new_point(points2, weights);
- Point<3> T = manifold.get_new_point(points3, weights);
+ Point<3> Q = manifold.get_new_point(make_array_view(points1),
+ make_array_view(weights));
+ Point<3> S = manifold.get_new_point(make_array_view(points2),
+ make_array_view(weights));
+ Point<3> T = manifold.get_new_point(make_array_view(points3),
+ make_array_view(weights));
Point<3> P5(0.707107, 0.707107, 0.0);
Point<3> P4(0.0, 0.0, 1.0);
// Quadrature (const std::vector< Point< dim > > &points, const std::vector< double > &weights);
return 0;
}
-
-
-
w[0] = 1.0-(double)i/((double)n_intermediates);
w[1] = 1.0 - w[0];
- Point<spacedim> ip = manifold.get_new_point(sp, w);
+ Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+ make_array_view(w));
Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
return 0;
}
-
std::vector<double> weights(2);
weights[0] = 0.1;
weights[1] = 0.9;
- Point<spacedim> p = cell->get_manifold().get_new_point(points, weights);
- Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+ Point<spacedim> p = cell->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
+ Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
deallog << "Distance between cell manifold and face manifold: "
<< (pref-p) << std::endl;
weights[0] = 0.55;
weights[1] = 0.45;
- p = cell->get_manifold().get_new_point(points, weights);
- pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+ p = cell->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
+ pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
deallog << "Distance between cell manifold and face manifold: "
<< (pref-p) << std::endl;
}
return 0;
}
-
std::vector<double> weights(2);
weights[0] = 0.1;
weights[1] = 0.9;
- Point<spacedim> p = cell->get_manifold().get_new_point(points, weights);
- Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+ Point<spacedim> p = cell->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
+ Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
deallog << "Distance between cell manifold and face manifold: "
<< (pref-p) << std::endl;
weights[0] = 0.55;
weights[1] = 0.45;
- p = cell->get_manifold().get_new_point(points, weights);
- pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+ p = cell->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
+ pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
deallog << "Distance between cell manifold and face manifold: "
<< (pref-p) << std::endl;
}
return 0;
}
-
std::vector<double> weights(2);
weights[0] = 0.1;
weights[1] = 0.9;
- Point<spacedim> p = cell->get_manifold().get_new_point(points, weights);
- Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+ Point<spacedim> p = cell->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
+ Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
deallog << "Distance between cell manifold and face manifold: "
<< (pref-p) << std::endl;
weights[0] = 0.55;
weights[1] = 0.45;
- p = cell->get_manifold().get_new_point(points, weights);
- pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+ p = cell->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
+ pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+ make_array_view(weights));
deallog << "Distance between cell manifold and face manifold: "
<< (pref-p) << std::endl;
}
return 0;
}
-
-