#include <deal.II/base/config.h>
+#include <deal.II/base/std_cxx11/array.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>
#include <deal.II/base/table.h>
{
/**
* @name Creating meshes for basic geometries
- * @{
*/
+ ///@{
/**
* Initialize the given triangulation with a hypercube (line in 1D, square
*
* @image html hyper_cubes.png
*
+ * If @p dim < @p spacedim, this will create a @p dim dimensional object
+ * in the first @p dim coordinate directions embedded into the @p spacedim
+ * dimensional space with the remaining entries set to zero. For example,
+ * a <tt>Triangulation@<2,3@></tt> will be a square in the xy plane with z=0.
+ *
* See also subdivided_hyper_cube() for a coarse mesh consisting of several
* cells. See hyper_rectangle(), if different lengths in different ordinate
* directions are required.
const bool colorize= false);
/**
- * \brief %Triangulation of a d-simplex with (d+1) vertices and mesh cells,
- * resp.
+ * \brief %Triangulation of a d-simplex with (d+1) vertices and mesh cells.
*
* The @p vertices argument contains a vector with all d+1 vertices of the
* simplex. They must be given in an order such that the vectors from the
* @image html simplex_2d.png
* @image html simplex_3d.png
*
+ * @param tria The Triangulation to create. It needs to be empty upon calling this
+ * function.
+ *
+ * @param vertices The dim+1 corners of the simplex.
+ *
+ * @note Implemented for <tt>Triangulation@<2,2@></tt>, <tt>Triangulation@<3,3@></tt>.
+ *
* @author Guido Kanschat
* @date 2015
*/
* cells. Thus, the number of cells filling the given volume is
* <tt>repetitions<sup>dim</sup></tt>.
*
- * If spacedim=dim+1 the same mesh as in the case spacedim=dim is created,
- * but the vertices have an additional coordinate =0. So, if dim=1 one
- * obtains line along the x axis in the xy plane, and if dim=3 one obtains a
- * square in lying in the xy plane in 3d space.
+ * If @p dim < @p spacedim, this will create a @p dim dimensional object
+ * in the first @p dim coordinate directions embedded into the @p spacedim
+ * dimensional space with the remaining entries set to zero. For example,
+ * a <tt>Triangulation@<2,3@></tt> will be a square in the xy plane with z=0.
*
* @note The triangulation needs to be void upon calling this function.
*/
* 2<sup>i</sup>. For instance, the center point (1,-1,1) yields a material
* id 5.
*
- * @note If spacedim>dim the same mesh as in the case spacedim=dim is
- * created, but the vertices have all additional coordinates equal to the
- * coordinates of p1, i.e., we generate a hyper_rectangle parallel to the
- * x-axis (1d), or to the xy-axis passing through p1. The additional
- * coordinates of the point p2 are ignored by this class in this case.
+ * If @p dim < @p spacedim, this will create a @p dim dimensional object
+ * in the first @p dim coordinate directions embedded into the @p spacedim
+ * dimensional space with the remaining entries set to zero. For example,
+ * a <tt>Triangulation@<2,3@></tt> will be a rectangle in the xy plane with z=0,
+ * defined by the two opposing corners @p p1 and @p p2.
*
* @note The triangulation needs to be void upon calling this function.
*/
template <int dim, int spacedim>
void hyper_rectangle (Triangulation<dim,spacedim> &tria,
- const Point<spacedim> &p1,
- const Point<spacedim> &p2,
+ const Point<dim> &p1,
+ const Point<dim> &p2,
const bool colorize = false);
/**
* Create a coordinate-parallel parallelepiped from the two diagonally
- * opposite corner points @p p1 and @p p2. In dimension @p i,
+ * opposite corner points @p p1 and @p p2. In direction @p i,
* <tt>repetitions[i]</tt> cells are generated.
*
* To get cells with an aspect ratio different from that of the domain, use
* be true. That means the boundary indicator is 0 on the left and 1 on the
* right. See step-15 for details.
*
- * @note The triangulation needs to be void upon calling this function.
+ * If @p dim < @p spacedim, this will create a @p dim dimensional object
+ * in the first @p dim coordinate directions embedded into the @p spacedim
+ * dimensional space with the remaining entries set to zero. For example,
+ * a <tt>Triangulation@<2,3@></tt> will be a rectangle in the xy plane with z=0,
+ * defined by the two opposing corners @p p1 and @p p2.
*
* @note For an example of the use of this function see the step-28 tutorial
* program.
+ *
+ * @param tria The Triangulation to create. It needs to be empty upon calling this
+ * function.
+ *
+ * @param repetitions A vector of dim positive values denoting the number of cells
+ * to generate in that direction.
+ *
+ * @param p1 First corner point.
+ *
+ * @param p2 Second corner opposite to @p p1.
+ *
+ * @param colorize Assign different boundary ids if set to true.
+ *
*/
template <int dim, int spacedim>
void
subdivided_hyper_rectangle (Triangulation<dim,spacedim> &tria,
const std::vector<unsigned int> &repetitions,
- const Point<spacedim> &p1,
- const Point<spacedim> &p2,
+ const Point<dim> &p1,
+ const Point<dim> &p2,
const bool colorize=false);
/**
*
* @image html cheese_2d.png
*
+ * If @p dim < @p spacedim, this will create a @p dim dimensional object
+ * in the first @p dim coordinate directions embedded into the @p spacedim
+ * dimensional space with the remaining entries set to zero.
+ *
+ * @param tria The Triangulation to create. It needs to be empty upon calling this
+ * function.
+ *
+ * @param holes Positive number of holes in each of the dim directions.
+
* @author Guido Kanschat
* @date 2015
*/
const Point<dim> (&corners) [dim],
const bool colorize = false);
+ /**
+ * A subdivided parallelepiped.
+ *
+ * @param tria The Triangulation to create. It needs to be empty upon calling this
+ * function.
+ *
+ * @param origin First corner of the parallelepiped.
+ *
+ * @param edges An array of @p dim tensors describing the length and direction of the edges
+ * from @p origin.
+ *
+ * @param subdivisions Number of subdivisions in each of the dim directions. Each
+ * entry must be positive. An empty vector is equivalent to one subdivision in
+ * each direction.
+ *
+ * @param colorize Assign different boundary ids if set to true.
+ *
+ * @note Implemented for all combinations of @p dim and @p spacedim.
+ */
+ template <int dim, int spacedim>
+ void
+ subdivided_parallelepiped (Triangulation<dim, spacedim> &tria,
+ const Point<spacedim> &origin,
+ const std_cxx11::array<Tensor<1,spacedim>,dim> &edges,
+ const std::vector<unsigned int> &subdivisions = std::vector<unsigned int>(),
+ const bool colorize = false);
/**
* Hypercube with a layer of hypercubes around it. The first two parameters
* This function is declared to exist for triangulations of all space
* dimensions, but throws an error if called in 1d.
*
+ * You should attach a SphericalManifold to the cells and faces for correct
+ * placement of vertices upon refinement and to be able to use higher order
+ * mappings.
+ *
* @note The triangulation needs to be void upon calling this function.
*/
template <int dim>
const double radius = 1.);
/**
- * Creates an hyper sphere.
- * This function is declared to exist for dim=1, spacedim =2, and
- * dim=2, spacedim =3.
+ * Creates a hyper sphere, i.e., a surface of a ball in @p spacedim
+ * dimensions.
+ * This function only exists for dim+1=spacedim in 2 and 3 space
+ * dimensions.
+ *
+ * You should attach a SphericalManifold to the cells and faces for correct
+ * placement of vertices upon refinement and to be able to use higher order
+ * mappings.
*
* The following pictures are generated with:
* @code
* triangulation.refine_global(3);
* @endcode
*
- * See the @ref manifold "documentation module on manifolds", for
- * more detail.
+ * See the @ref manifold "documentation module on manifolds" for
+ * more details.
*
* @image html sphere.png
* @image html sphere_section.png
const double R,
const double r);
- /*
- * @}
- */
+ ///@}
/**
* @name Creating meshes from other meshes
- * @{
*/
+ ///@{
/**
* Given the two triangulations specified as the first two arguments, create
void flatten_triangulation(const Triangulation<dim,spacedim1> &in_tria,
Triangulation<dim,spacedim2> &out_tria);
- /*
- * @}
- */
+ ///@}
/**
- * @name Creating lower-dimensional meshes from parts of higher-dimensional
- * meshes
- * @{
+ * @name Creating lower-dimensional meshes from parts of higher-dimensional meshes
*/
-
+ ///@{
#ifdef _MSC_VER
// Microsoft's VC++ has a bug where it doesn't want to recognize that
const std::set<types::boundary_id> &boundary_ids
= std::set<types::boundary_id>());
- /*
- * @}
- */
+ ///@}
+
/**
* @name Exceptions
- * @{
*/
+ ///@{
+
/**
* Exception
<< "The vector of repetitions must have "
<< arg1 <<" elements.");
- /*
- * @}
- */
+ ///@}
}
template <int dim, int spacedim>
void
hyper_rectangle (Triangulation<dim,spacedim> &tria,
- const Point<spacedim> &p_1,
- const Point<spacedim> &p_2,
+ const Point<dim> &p_1,
+ const Point<dim> &p_2,
const bool colorize)
{
- // First, normalize input such that
- // p1 is lower in all coordinate directions.
- Point<spacedim> p1(p_1);
- Point<spacedim> p2(p_2);
-
- for (unsigned int i=0; i<spacedim; ++i)
- if (p1(i) > p2(i))
- std::swap (p1(i), p2(i));
+ // First, extend dimensions from dim to spacedim and
+ // normalize such that p1 is lower in all coordinate
+ // directions. Additional entries will be 0.
+ Point<spacedim> p1, p2;
+ for (unsigned int i=0; i<dim; ++i)
+ {
+ p1(i) = std::min(p_1(i), p_2(i));
+ p2(i) = std::max(p_1(i), p_2(i));
+ }
std::vector<Point<spacedim> > vertices (GeometryInfo<dim>::vertices_per_cell);
switch (dim)
Assert (left < right,
ExcMessage ("Invalid left-to-right bounds of hypercube"));
- Point<spacedim> p1;
- Point<spacedim> p2;
-
- p1(spacedim-1) = 0;
- p2(spacedim-1) = 0;
-
+ Point<dim> p1, p2;
for (unsigned int i=0; i<dim; ++i)
{
p1(i) = left;
const Point<2> (&corners)[2],
const bool colorize)
{
- std::vector<Point<2> > vertices (GeometryInfo<2>::vertices_per_cell);
-
- vertices[1] = corners[0];
- vertices[2] = corners[1];
- vertices[3] = vertices[1] + vertices[2];
- // Prepare cell data
- std::vector<CellData<2> > cells (1);
- for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
- cells[0].vertices[i] = i;
- cells[0].material_id = 0;
-
- tria.create_triangulation (vertices, cells, SubCellData());
-
- // Assign boundary indicators
- if (colorize)
- colorize_hyper_rectangle (tria);
+ Point<2> origin;
+ std_cxx11::array<Tensor<1,2>,2> edges;
+ edges[0] = corners[0];
+ edges[1] = corners[1];
+ std::vector<unsigned int> subdivisions;
+ subdivided_parallelepiped(tria, origin, edges, subdivisions, colorize);
}
colorize);
}
- // Parallelepiped implementation in 1d, 2d, and 3d. @note The
- // implementation in 1d is similar to hyper_rectangle(), and in 2d is
- // similar to parallelogram().
- //
- // The GridReordering::reorder_grid is made use of towards the end of
- // this function. Thus the triangulation is explicitly constructed for
- // all dim here since it is slightly different in that respect
- // (cf. hyper_rectangle(), parallelogram()).
template<int dim>
void
subdivided_parallelepiped (Triangulation<dim> &tria,
const Point<dim> (&corners) [dim],
const bool colorize)
{
- // Zero n_subdivisions is the origin only, which makes no sense
+ Point<dim> origin;
+ std::vector<unsigned int> subdivisions;
+ std_cxx11::array<Tensor<1,dim>,dim> edges;
+ for (unsigned int i=0; i<dim; ++i)
+ {
+ subdivisions.push_back(n_subdivisions[i]);
+ edges[i] = corners[i];
+ }
+
+ subdivided_parallelepiped (tria, origin, edges, subdivisions, colorize);
+ }
+
+ // Parallelepiped implementation in 1d, 2d, and 3d. @note The
+ // implementation in 1d is similar to hyper_rectangle(), and in 2d is
+ // similar to parallelogram().
+ //
+ // The GridReordering::reorder_grid is made use of towards the end of
+ // this function. Thus the triangulation is explicitly constructed for
+ // all dim here since it is slightly different in that respect
+ // (cf. hyper_rectangle(), parallelogram()).
+ template <int dim, int spacedim>
+ void
+ subdivided_parallelepiped (Triangulation<dim, spacedim> &tria,
+ const Point<spacedim> &origin,
+ const std_cxx11::array<Tensor<1,spacedim>,dim> &edges,
+ const std::vector<unsigned int> &subdivisions,
+ const bool colorize)
+ {
+ if (subdivisions.size()==0)
+ {
+ std::vector<unsigned int> new_subdivisions(dim, 1);
+ subdivided_parallelepiped(tria, origin, edges, new_subdivisions, colorize);
+ return;
+ }
+
+ Assert(subdivisions.size()==dim, ExcMessage(""));
+
+ // check subdivisions
for (unsigned int i=0; i<dim; ++i)
- Assert (n_subdivisions[i]>0, ExcInvalidRepetitions(n_subdivisions[i]));
+ {
+ Assert (subdivisions[i]>0, ExcInvalidRepetitions(subdivisions[i]));
+ Assert (edges[i].norm()>0, ExcMessage("Edges in subdivided_parallelepiped() must not be degenerate."));
+ }
// Check corners do not overlap (unique)
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=i+1; j<dim; ++j)
- Assert ((corners[i]!=corners[j]),
- ExcMessage ("Invalid distance between corner points of parallelepiped."));
+ Assert ((edges[i]!=edges[j]),
+ ExcMessage ("Degenerate edges of subdivided_parallelepiped encountered."));
// Create a list of points
- Point<dim> delta[dim];
-
- for (unsigned int i=0; i<dim; ++i)
- delta[i] = corners[i]/n_subdivisions[i];
- std::vector<Point<dim> > points;
+ std::vector<Point<spacedim> > points;
switch (dim)
{
case 1:
- for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
- points.push_back (double(x)*delta[0]);
+ for (unsigned int x=0; x<=subdivisions[0]; ++x)
+ points.push_back (origin + edges[0]/subdivisions[0]*x);
break;
case 2:
- for (unsigned int y=0; y<=n_subdivisions[1]; ++y)
- for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
- points.push_back (double(x)*delta[0] + double(y)*delta[1]);
+ for (unsigned int y=0; y<=subdivisions[1]; ++y)
+ for (unsigned int x=0; x<=subdivisions[0]; ++x)
+ points.push_back (origin
+ + edges[0]/subdivisions[0]*x
+ + edges[1]/subdivisions[1]*y);
break;
case 3:
- for (unsigned int z=0; z<=n_subdivisions[2]; ++z)
- for (unsigned int y=0; y<=n_subdivisions[1]; ++y)
- for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
- points.push_back (double(x)*delta[0] + double(y)*delta[1] + double(z)*delta[2]);
+ for (unsigned int z=0; z<=subdivisions[2]; ++z)
+ for (unsigned int y=0; y<=subdivisions[1]; ++y)
+ for (unsigned int x=0; x<=subdivisions[0]; ++x)
+ points.push_back (
+ origin
+ + edges[0]/subdivisions[0]*x
+ + edges[1]/subdivisions[1]*y
+ + edges[2]/subdivisions[2]*z);
break;
default:
// Prepare cell data
unsigned int n_cells = 1;
for (unsigned int i=0; i<dim; ++i)
- n_cells *= n_subdivisions[i];
+ n_cells *= subdivisions[i];
std::vector<CellData<dim> > cells (n_cells);
// Create fixed ordering of
switch (dim)
{
case 1:
- for (unsigned int x=0; x<n_subdivisions[0]; ++x)
+ for (unsigned int x=0; x<subdivisions[0]; ++x)
{
cells[x].vertices[0] = x;
cells[x].vertices[1] = x+1;
case 2:
{
// Shorthand
- const unsigned int n_dy = n_subdivisions[1];
- const unsigned int n_dx = n_subdivisions[0];
+ const unsigned int n_dy = subdivisions[1];
+ const unsigned int n_dx = subdivisions[0];
for (unsigned int y=0; y<n_dy; ++y)
for (unsigned int x=0; x<n_dx; ++x)
{
- const unsigned int c = y*n_dx + x;
+ const unsigned int c = y*n_dx + x;
cells[c].vertices[0] = y*(n_dx+1) + x;
cells[c].vertices[1] = y*(n_dx+1) + x+1;
cells[c].vertices[2] = (y+1)*(n_dx+1) + x;
case 3:
{
// Shorthand
- const unsigned int n_dz = n_subdivisions[2];
- const unsigned int n_dy = n_subdivisions[1];
- const unsigned int n_dx = n_subdivisions[0];
+ const unsigned int n_dz = subdivisions[2];
+ const unsigned int n_dy = subdivisions[1];
+ const unsigned int n_dx = subdivisions[0];
for (unsigned int z=0; z<n_dz; ++z)
for (unsigned int y=0; y<n_dy; ++y)
for (unsigned int x=0; x<n_dx; ++x)
{
- const unsigned int c = z*n_dy*n_dx + y*n_dx + x;
+ const unsigned int c = z*n_dy*n_dx + y*n_dx + x;
cells[c].vertices[0] = z*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x;
cells[c].vertices[1] = z*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x+1;
Assert (left < right,
ExcMessage ("Invalid left-to-right bounds of hypercube"));
- Point<spacedim> p0,p1;
+ Point<dim> p0, p1;
for (unsigned int i=0; i<dim; ++i)
{
p0[i] = left;
subdivided_hyper_rectangle (
Triangulation<dim, spacedim> &tria,
const std::vector<unsigned int> &repetitions,
- const Point<spacedim> &p_1,
- const Point<spacedim> &p_2,
+ const Point<dim> &p_1,
+ const Point<dim> &p_2,
const bool colorize)
{
- // contributed by Joerg R. Weimar
- // (j.weimar@jweimar.de) 2003
Assert(repetitions.size() == dim,
ExcInvalidRepetitionsDimension(dim));
- // First, normalize input such that
- // p1 is lower in all coordinate
- // directions.
- Point<spacedim> p1(p_1);
- Point<spacedim> p2(p_2);
+ // First, extend dimensions from dim to spacedim and
+ // normalize such that p1 is lower in all coordinate
+ // directions. Additional entries will be 0.
+ Point<spacedim> p1, p2;
for (unsigned int i=0; i<dim; ++i)
- if (p1(i) > p2(i))
- std::swap (p1(i), p2(i));
+ {
+ p1(i) = std::min(p_1(i), p_2(i));
+ p2(i) = std::max(p_1(i), p_2(i));
+ }
- // then check that all repetitions
- // are >= 1, and calculate deltas
- // convert repetitions from double
- // to int by taking the ceiling.
+ // calculate deltas and validate input
std::vector<Point<spacedim> > delta(dim);
-
for (unsigned int i=0; i<dim; ++i)
{
Assert (repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
delta[i][i] = (p2[i]-p1[i])/repetitions[i];
+ Assert(delta[i][i]>0.0,
+ ExcMessage("The first dim entries of coordinates of p1 and p2 need to be different."));
}
- // then generate the necessary
- // points
+ // then generate the points
std::vector<Point<spacedim> > points;
switch (dim)
{
}
// next create the cells
- // Prepare cell data
std::vector<CellData<dim> > cells;
switch (dim)
{
const Point<dim> &p_2,
const bool colorize)
{
- // contributed by Joerg R. Weimar
- // (j.weimar@jweimar.de) 2003
- // modified by Yaqi Wang 2006
Assert(step_sz.size() == dim,
ExcInvalidRepetitionsDimension(dim));
-
// First, normalize input such that
// p1 is lower in all coordinate
- // directions.
-
- // and check the consistency of
+ // directions and check the consistency of
// step sizes, i.e. that they all
// add up to the sizes specified by
// p_1 and p_2
const Table<1,types::material_id> &material_id,
const bool colorize)
{
- // contributed by Yaqi Wang 2006
Assert(spacing.size() == 1,
ExcInvalidRepetitionsDimension(1));
const Table<2,types::material_id> &material_id,
const bool colorize)
{
- // contributed by Yaqi Wang 2006
Assert(spacing.size() == 2,
ExcInvalidRepetitionsDimension(2));
const Table<3,types::material_id> &material_id,
const bool colorize)
{
- // contributed by Yaqi Wang 2006
const unsigned int dim = 3;
Assert(spacing.size() == dim,
template void
hyper_rectangle<deal_II_dimension, deal_II_space_dimension> (
Triangulation<deal_II_dimension, deal_II_space_dimension> &,
- const Point<deal_II_space_dimension>&, const Point<deal_II_space_dimension>&,
+ const Point<deal_II_dimension>&, const Point<deal_II_dimension>&,
const bool);
template void
subdivided_hyper_rectangle<deal_II_dimension, deal_II_space_dimension>
(Triangulation<deal_II_dimension, deal_II_space_dimension> &,
const std::vector<unsigned int>&,
- const Point<deal_II_space_dimension>&,
- const Point<deal_II_space_dimension>&,
+ const Point<deal_II_dimension>&,
+ const Point<deal_II_dimension>&,
const bool);
+ template void
+ subdivided_parallelepiped<deal_II_dimension, deal_II_space_dimension>
+ (Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+ const Point<deal_II_space_dimension> &,
+ const std_cxx11::array<Tensor<1,deal_II_space_dimension>,deal_II_dimension> &,
+ const std::vector<unsigned int> &,
+ const bool colorize);
+
template void
hyper_cross<deal_II_dimension, deal_II_space_dimension>
(Triangulation<deal_II_dimension, deal_II_space_dimension> &,
const Point<deal_II_dimension> (&) [deal_II_dimension],
const bool);
#endif
+
\}
}
template<int dim, int spacedim>
void test(std::ostream &out)
{
- Point<spacedim> p1;
+ Point<dim> p1;
p1[0] = 2.;
if (dim>1) p1[1] = -1.;
- if (dim>2) p1[2] = 0.;
- Point<spacedim> p2;
+ Point<dim> p2;
p2[0] = 3.;
if (dim>1) p2[1] = 2.;
- if (dim>2) p2[2] = 4.;
- Point<spacedim> p3;
+ Point<dim> p3;
p3[0] = 2.;
if (dim>1) p3[1] = 1.;
- if (dim>2) p3[2] = 4.;
GridOut go;
GridOut::OutputFormat format = GridOut::msh;
}
}
-
int main()
{
std::ofstream logfile("output");
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test subdivided_parallelepiped
+
+#include "../tests.h"
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/logstream.h>
+#include <cmath>
+#include <cstdlib>
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+template <int dim>
+Point<dim> point(double x=0, double y=0, double z=0)
+{
+ Point<dim> p;
+ if (dim>0) p[0] = x;
+ if (dim>1) p[1] = y;
+ if (dim>2) p[2] = z;
+ return p;
+}
+
+
+// The simplest test case is to create a parallelepiped grid with a
+// number of subdivisions and output the result.
+template<int dim, int spacedim>
+void check (bool subdivide)
+{
+ deallog << "dim " << dim << " spacedim " << spacedim << std::endl;
+
+ Point<spacedim> origin;
+ for (unsigned int d=0; d<spacedim; ++d)
+ origin[d] = 0.1+d*1.0;
+
+ std_cxx11::array<Tensor<1,spacedim>,dim> edges;
+ switch (dim)
+ {
+ case 1:
+ edges[0] = point<spacedim>(0.5, 0.02, 0.03);
+ break;
+
+ case 2:
+ edges[0] = point<spacedim>(0.15, 0.50, 0.03);
+ edges[1] = point<spacedim>(0.70, 0.25, -0.01);
+ break;
+
+ case 3:
+ edges[0] = point<spacedim>(0.10, 0.50, 0.70);
+ edges[1] = point<spacedim>(1.50, 0.25, 0.70);
+ edges[2] = point<spacedim>(0.25, 0.50, 0.3);
+ break;
+
+ default:
+ Assert (false, ExcInternalError ());
+ }
+
+
+ std::vector<unsigned int> subdivisions;
+ if (subdivide)
+ {
+ subdivisions.resize(dim);
+ for (unsigned int d=0; d<dim; ++d)
+ subdivisions[d]=d+1;
+ }
+
+
+ bool colorize = false;
+
+ Triangulation<dim,spacedim> triangulation;
+ GridGenerator::subdivided_parallelepiped (triangulation, origin, edges,
+ subdivisions, colorize);
+
+ GridOut grid_out;
+
+ grid_out.write_gnuplot (triangulation, deallog.get_file_stream());
+}
+
+int main ()
+{
+ initlog();
+
+ check<1,1> (true);
+ check<1,2> (true);
+ check<1,3> (true);
+
+ check<2,2> (true);
+ check<2,3> (true);
+
+ check<3,3> (true);
+ check<3,3> (false);
+
+
+ deallog<< "OK" << std::endl;
+}
--- /dev/null
+
+DEAL::dim 1 spacedim 1
+0.100000 0 0
+0.600000 0 0
+
+
+DEAL::dim 1 spacedim 2
+0.100000 1.10000 0 0
+0.600000 1.12000 0 0
+
+
+DEAL::dim 1 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.600000 1.12000 2.13000 0 0
+
+
+DEAL::dim 2 spacedim 2
+0.100000 1.10000 0 0
+0.250000 1.60000 0 0
+0.600000 1.72500 0 0
+0.450000 1.22500 0 0
+0.100000 1.10000 0 0
+
+
+0.450000 1.22500 0 0
+0.600000 1.72500 0 0
+0.950000 1.85000 0 0
+0.800000 1.35000 0 0
+0.450000 1.22500 0 0
+
+
+DEAL::dim 2 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.250000 1.60000 2.13000 0 0
+0.600000 1.72500 2.12500 0 0
+0.450000 1.22500 2.09500 0 0
+0.100000 1.10000 2.10000 0 0
+
+
+0.450000 1.22500 2.09500 0 0
+0.600000 1.72500 2.12500 0 0
+0.950000 1.85000 2.12000 0 0
+0.800000 1.35000 2.09000 0 0
+0.450000 1.22500 2.09500 0 0
+
+
+DEAL::dim 3 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.200000 1.60000 2.80000 0 0
+0.283333 1.76667 2.90000 0 0
+0.183333 1.26667 2.20000 0 0
+0.100000 1.10000 2.10000 0 0
+
+0.850000 1.22500 2.45000 0 0
+0.950000 1.72500 3.15000 0 0
+1.03333 1.89167 3.25000 0 0
+0.933333 1.39167 2.55000 0 0
+0.850000 1.22500 2.45000 0 0
+
+0.100000 1.10000 2.10000 0 0
+0.850000 1.22500 2.45000 0 0
+
+0.200000 1.60000 2.80000 0 0
+0.950000 1.72500 3.15000 0 0
+
+0.283333 1.76667 2.90000 0 0
+1.03333 1.89167 3.25000 0 0
+
+0.183333 1.26667 2.20000 0 0
+0.933333 1.39167 2.55000 0 0
+
+0.850000 1.22500 2.45000 0 0
+0.950000 1.72500 3.15000 0 0
+1.03333 1.89167 3.25000 0 0
+0.933333 1.39167 2.55000 0 0
+0.850000 1.22500 2.45000 0 0
+
+1.60000 1.35000 2.80000 0 0
+1.70000 1.85000 3.50000 0 0
+1.78333 2.01667 3.60000 0 0
+1.68333 1.51667 2.90000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.850000 1.22500 2.45000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.950000 1.72500 3.15000 0 0
+1.70000 1.85000 3.50000 0 0
+
+1.03333 1.89167 3.25000 0 0
+1.78333 2.01667 3.60000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.68333 1.51667 2.90000 0 0
+
+0.183333 1.26667 2.20000 0 0
+0.283333 1.76667 2.90000 0 0
+0.366667 1.93333 3.00000 0 0
+0.266667 1.43333 2.30000 0 0
+0.183333 1.26667 2.20000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.03333 1.89167 3.25000 0 0
+1.11667 2.05833 3.35000 0 0
+1.01667 1.55833 2.65000 0 0
+0.933333 1.39167 2.55000 0 0
+
+0.183333 1.26667 2.20000 0 0
+0.933333 1.39167 2.55000 0 0
+
+0.283333 1.76667 2.90000 0 0
+1.03333 1.89167 3.25000 0 0
+
+0.366667 1.93333 3.00000 0 0
+1.11667 2.05833 3.35000 0 0
+
+0.266667 1.43333 2.30000 0 0
+1.01667 1.55833 2.65000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.03333 1.89167 3.25000 0 0
+1.11667 2.05833 3.35000 0 0
+1.01667 1.55833 2.65000 0 0
+0.933333 1.39167 2.55000 0 0
+
+1.68333 1.51667 2.90000 0 0
+1.78333 2.01667 3.60000 0 0
+1.86667 2.18333 3.70000 0 0
+1.76667 1.68333 3.00000 0 0
+1.68333 1.51667 2.90000 0 0
+
+0.933333 1.39167 2.55000 0 0
+1.68333 1.51667 2.90000 0 0
+
+1.03333 1.89167 3.25000 0 0
+1.78333 2.01667 3.60000 0 0
+
+1.11667 2.05833 3.35000 0 0
+1.86667 2.18333 3.70000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.76667 1.68333 3.00000 0 0
+
+0.266667 1.43333 2.30000 0 0
+0.366667 1.93333 3.00000 0 0
+0.450000 2.10000 3.10000 0 0
+0.350000 1.60000 2.40000 0 0
+0.266667 1.43333 2.30000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.11667 2.05833 3.35000 0 0
+1.20000 2.22500 3.45000 0 0
+1.10000 1.72500 2.75000 0 0
+1.01667 1.55833 2.65000 0 0
+
+0.266667 1.43333 2.30000 0 0
+1.01667 1.55833 2.65000 0 0
+
+0.366667 1.93333 3.00000 0 0
+1.11667 2.05833 3.35000 0 0
+
+0.450000 2.10000 3.10000 0 0
+1.20000 2.22500 3.45000 0 0
+
+0.350000 1.60000 2.40000 0 0
+1.10000 1.72500 2.75000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.11667 2.05833 3.35000 0 0
+1.20000 2.22500 3.45000 0 0
+1.10000 1.72500 2.75000 0 0
+1.01667 1.55833 2.65000 0 0
+
+1.76667 1.68333 3.00000 0 0
+1.86667 2.18333 3.70000 0 0
+1.95000 2.35000 3.80000 0 0
+1.85000 1.85000 3.10000 0 0
+1.76667 1.68333 3.00000 0 0
+
+1.01667 1.55833 2.65000 0 0
+1.76667 1.68333 3.00000 0 0
+
+1.11667 2.05833 3.35000 0 0
+1.86667 2.18333 3.70000 0 0
+
+1.20000 2.22500 3.45000 0 0
+1.95000 2.35000 3.80000 0 0
+
+1.10000 1.72500 2.75000 0 0
+1.85000 1.85000 3.10000 0 0
+
+DEAL::dim 3 spacedim 3
+0.100000 1.10000 2.10000 0 0
+0.200000 1.60000 2.80000 0 0
+0.450000 2.10000 3.10000 0 0
+0.350000 1.60000 2.40000 0 0
+0.100000 1.10000 2.10000 0 0
+
+1.60000 1.35000 2.80000 0 0
+1.70000 1.85000 3.50000 0 0
+1.95000 2.35000 3.80000 0 0
+1.85000 1.85000 3.10000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.100000 1.10000 2.10000 0 0
+1.60000 1.35000 2.80000 0 0
+
+0.200000 1.60000 2.80000 0 0
+1.70000 1.85000 3.50000 0 0
+
+0.450000 2.10000 3.10000 0 0
+1.95000 2.35000 3.80000 0 0
+
+0.350000 1.60000 2.40000 0 0
+1.85000 1.85000 3.10000 0 0
+
+DEAL::OK
projects / dealii.git / commitdiff