// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* set, parts of the boundary receive different boundary numbers,
* allowing them to be distinguished by application programs. See the
* documentation of the functions for details.
- *
+ *
* Additionally this class provides a function
* (@p laplace_transformation) that smoothly transforms a grid
* according to given new boundary points. This can be used to
* adds 2<sup>i</sup>. For
* instance, the center point
* (1,-1,1) yields a material id 5.
- *
+ *
* @note The triangulation needs to be
* void upon calling this
* function.
* dimension @p i,
* <tt>repetitions[i]</tt> cells are
* generated.
- *
+ *
* To get cells with an aspect
* ratio different from that of
* the domain, use different
* of integers denoting the
* number of subdivisions in each
* coordinate direction.
- *
+ *
* If the @p colorize flag is
* set, the
* @p boundary_indicators of the
const Point<dim> &p,
const Table<dim,unsigned char> &material_id,
const bool colorize=false);
-
+
/**
* A parallelogram. The first
* corner point is the
parallelogram(Triangulation<dim>& tria,
const Tensor<2,dim>& corners,
const bool colorize=false);
-
-
+
+
/**
* Hypercube with a layer of
* hypercubes around it. The
const double right= 1.,
const double thickness = 1.,
const bool colorize = false);
-
+
/**
* Initialize the given
* triangulation with a
* @note The triangulation needs to be
* void upon calling this
* function.
- */
+ */
template <int dim>
static void hyper_ball (Triangulation<dim> &tria,
const Point<dim> ¢er = Point<dim>(),
* perpendicular to the
* <i>x</i>-axis.
*
- * The boundary indicators for the final
+ * The boundary indicators for the final
* triangulation are 0 for the curved boundary and
- * 1 for the cut plane.
+ * 1 for the cut plane.
*
* The appropriate
- * boundary class is
+ * boundary class is
* HalfHyperBallBoundary, or HyperBallBoundary.
- *
+ *
* @note The triangulation needs to be
* void upon calling this
* function.
static void cylinder (Triangulation<dim> &tria,
const double radius = 1.,
const double half_length = 1.);
-
+
+ /**
+ * Create a cutted cone around
+ * the x-axis. The cone extends
+ * from <tt>x=-half_length</tt>
+ * to <tt>x=half_length</tt> and
+ * its projection into the @p
+ * yz-plane is a circle of radius
+ * @p radius_0 at
+ * <tt>x=-half_length</tt> and a
+ * circle of radius @p radius_1
+ * at <tt>x=+half_length</tt>.
+ * In between the radius is
+ * linearly decreasing.
+ *
+ * In two dimensions, the cone is
+ * a trapezoid from
+ * <tt>x=-half_length</tt> to
+ * <tt>x=+half_length</tt> and
+ * from <tt>y=-radius</tt> to
+ * <tt>y=radius_0</tt> at
+ * <tt>x=-half_length</tt> and
+ * from <tt>y=-radius_1</tt> to
+ * <tt>y=radius_1</tt> at
+ * <tt>x=+half_length</tt>. In
+ * between the range of
+ * <tt>y</tt> is linearly
+ * decreasing.
+ *
+ * The boundaries are colored
+ * according to the following
+ * scheme: 0 for the hull of the
+ * cone, 1 for the left hand
+ * face and 2 for the right hand
+ * face.
+ *
+ * @note The triangulation needs to be
+ * void upon calling this
+ * function.
+ *
+ * @author Markus Bürg, 2009
+ */
+ template <int dim>
+ static void
+ truncated_cone (Triangulation<dim> &tria,
+ const double radius_0 = 1.0,
+ const double radius_1 = 0.5,
+ const double half_length = 1.0);
+
/**
* Initialize the given
* triangulation with a hyper-L
static void hyper_L (Triangulation<dim> &tria,
const double left = -1.,
const double right= 1.);
-
+
/**
* Initialize the given
* Triangulation with a hypercube
* <i>z</i>-direction, such that
* a plane cuts the lower half of
* a rectangle in two.
-
+
* This function is declared to
* exist for triangulations of
* all space dimensions, but
const double left = 0.,
const double right= 1.,
const bool colorize = false);
-
+
/**
* Produce a hyper-shell,
* the region between two
const double outer_radius,
const unsigned int n_cells = 0,
bool colorize = false);
-
+
/**
* Produce a half hyper-shell,
* i.e. the space between two
const double inner_radius,
const double outer_radius,
const unsigned int n_cells = 0);
-
+
/**
* Produce a domain that is the space
* between two cylinders in 3d, with
const unsigned int n_radial_cells = 0,
const unsigned int n_axial_cells = 0);
- /**
+ /**
* This class produces a square
* on the <i>xy</i>-plane with a
* circular hole in the middle,
* (only in 3d).
*
* @image html cubes_hole.png
- *
+ *
* It is implemented in 2d and
* 3d, and takes the following
* arguments:
- *
+ *
* @arg @p inner_radius: size of the
- * internal hole
+ * internal hole
* @arg @p outer_radius: size of the
* biggest enclosed cylinder
* @arg @p L: extension on the @p z-direction
* @arg @p repetitions: number of subdivisions
* along the @p z-direction
- * @arg @p colorize: wether to assign different
+ * @arg @p colorize: wether to assign different
* boundary indicators to different faces.
- * The colors are given in lexicographic
- * ordering for the flat faces (0 to 3 in 2d,
- * 0 to 5 in 3d) plus the curved hole
+ * The colors are given in lexicographic
+ * ordering for the flat faces (0 to 3 in 2d,
+ * 0 to 5 in 3d) plus the curved hole
* (4 in 2d, and 6 in 3d).
- * If @p colorize is set to false, then flat faces
+ * If @p colorize is set to false, then flat faces
* get the number 0 and the hole gets number 1.
*/
template<int dim>
- static void hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
+ static void hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
const double inner_radius = .25,
const double outer_radius = .5,
const double L = .5,
const unsigned int n_rotations,
const double R,
const double r);
-
+
/**
* This function transformes the
* @p Triangulation @p tria
*/
DeclException1 (ExcInvalidRepetitionsDimension,
int,
- << "The vector of repetitions must have "
+ << "The vector of repetitions must have "
<< arg1 <<" elements.");
private:
const Point<dim>& center,
const double inner_radius,
const double outer_radius);
-
+
/**
* Solve the Laplace equation for
* @p laplace_transformation
Point<3>( 0, 0,-1),
Point<3>( 0, 0, 1)
};
-
+
#endif
}
// 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));
-
+
std::vector<Point<spacedim> > vertices (GeometryInfo<dim>::vertices_per_cell);
switch (dim)
{
case 2:
vertices[0] = vertices[1] = p1;
vertices[2] = vertices[3] = p2;
-
+
vertices[1](0) = p2(0);
vertices[2](0) = p1(0);
break;
case 3:
vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
-
+
vertices[1](0) = p2(0);
vertices[2](1) = p2(1);
vertices[3](0) = p2(0);
vertices[3](1) = p2(1);
-
+
vertices[4](0) = p1(0);
vertices[4](1) = p1(1);
vertices[5](1) = p1(1);
vertices[6](0) = p1(0);
-
+
break;
default:
Assert (false, ExcNotImplemented ());
{
Assert (left < right,
ExcMessage ("Invalid left and right bounds of hypercube"));
-
+
Point<spacedim> p1;
Point<spacedim> p2;
Assert (n_cells>4, ExcMessage("More than 4 cells are needed to create a moebius grid."));
Assert (r>0 && R>0, ExcMessage("Outer and inner radius must be positive."));
Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
-
-
+
+
std::vector<Point<dim> > vertices (4*n_cells);
double beta_step=n_rotations*numbers::PI/2.0/n_cells;
double alpha_step=2.0*numbers::PI/n_cells;
-
+
for (unsigned int i=0; i<n_cells; ++i)
for (unsigned int j=0; j<4; ++j)
{
}
unsigned int offset=0;
-
+
std::vector<CellData<dim> > cells (n_cells);
for (unsigned int i=0; i<n_cells; ++i)
{
offset+=4;
cells[i].material_id=0;
}
-
+
// now correct the last four vertices
cells[n_cells-1].vertices[4]=(0+n_rotations)%4;
cells[n_cells-1].vertices[5]=(3+n_rotations)%4;
const bool colorize)
{
Assert (dim==2, ExcNotImplemented());
-
+
std::vector<Point<dim> > vertices (GeometryInfo<dim>::vertices_per_cell);
vertices[1] = corners[0];
Assert (repetitions >= 1, ExcInvalidRepetitions(repetitions));
Assert (left < right,
ExcMessage ("Invalid left and right bounds of hypercube"));
-
+
// first generate the necessary
// points
const double delta = (right-left)/repetitions;
cells[x].material_id = 0;
}
break;
-
+
case 2:
cells.resize (repetitions*repetitions);
for (unsigned int y=0; y<repetitions; ++y)
Assert (false, ExcNotImplemented());
}
- tria.create_triangulation (points, cells, SubCellData());
+ tria.create_triangulation (points, cells, SubCellData());
}
{
// contributed by Joerg R. Weimar
// (j.weimar@jweimar.de) 2003
- Assert(repetitions.size() == dim,
+ Assert(repetitions.size() == dim,
ExcInvalidRepetitionsDimension(dim));
// First, normalize input such that
// p1 is lower in all coordinate
// directions.
Point<dim> p1(p_1);
Point<dim> p2(p_2);
-
+
for (unsigned int i=0;i<dim;++i)
if (p1(i) > p2(i))
std::swap (p1(i), p2(i));
// convert repetitions from double
// to int by taking the ceiling.
Point<dim> delta;
-
+
for (unsigned int i=0; i<dim; ++i)
{
Assert (repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
-
+
delta[i] = (p2[i]-p1[i])/repetitions[i];
}
-
+
// then generate the necessary
// points
std::vector<Point<dim> > points;
}
break;
}
-
+
case 2:
{
cells.resize (repetitions[1]*repetitions[0]);
}
break;
}
-
+
case 3:
{
const unsigned int n_x = (repetitions[0]+1);
const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1);
-
+
cells.resize (repetitions[2]*repetitions[1]*repetitions[0]);
for (unsigned int z=0; z<repetitions[2]; ++z)
for (unsigned int y=0; y<repetitions[1]; ++y)
cells[c].material_id = 0;
}
break;
-
+
}
default:
Assert (false, ExcNotImplemented());
}
- tria.create_triangulation (points, cells, SubCellData());
+ tria.create_triangulation (points, cells, SubCellData());
if (colorize)
{
// faces of all cells and set
// boundary indicator to the
// correct value if it was 0.
-
+
// use a large epsilon to
// compare numbers to avoid
// roundoff problems.
= 0.01 * *std::min_element (&delta[0], &delta[0]+dim);
Assert (epsilon > 0,
ExcMessage ("The distance between corner points must be positive."))
-
+
// actual code is external since
// 1-D is different from 2/3D.
colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon);
// contributed by Joerg R. Weimar
// (j.weimar@jweimar.de) 2003
// modified by Yaqi Wang 2006
- Assert(step_sz.size() == dim,
+ Assert(step_sz.size() == dim,
ExcInvalidRepetitionsDimension(dim));
}
double x = 0;
- for (unsigned int j=0; j<step_sizes.at(i).size(); j++)
+ for (unsigned int j=0; j<step_sizes.at(i).size(); j++)
x += step_sizes[i][j];
Assert(std::fabs(x - (p2(i)-p1(i))) <= 1e-12*std::fabs(x),
ExcInvalidRepetitions (i) );
}
-
+
// then generate the necessary
// points
std::vector<Point<dim> > points;
{
double x=0;
for (unsigned int i=0; ; ++i)
- {
+ {
points.push_back (Point<dim> (p1[0]+x));
// form partial sums. in
// and exit early instead
if (i == step_sizes[0].size())
break;
-
+
x += step_sizes[0][i];
}
break;
}
-
+
case 2:
{
double y=0;
p1[1]+y));
if (i == step_sizes[0].size())
break;
-
+
x += step_sizes[0][i];
}
-
+
if (j == step_sizes[1].size())
break;
-
+
y += step_sizes[1][j];
}
break;
p1[2]+z));
if (i == step_sizes[0].size())
break;
-
+
x += step_sizes[0][i];
}
if (j == step_sizes[1].size())
break;
-
+
y += step_sizes[1][j];
}
-
+
if (k == step_sizes[2].size())
break;
-
+
z += step_sizes[2][k];
}
break;
}
-
+
default:
Assert (false, ExcNotImplemented());
}
}
break;
}
-
+
case 2:
{
cells.resize (step_sizes[1].size()*step_sizes[0].size());
}
break;
}
-
+
case 3:
{
const unsigned int n_x = (step_sizes[0].size()+1);
const unsigned int n_xy = (step_sizes[0].size()+1)*(step_sizes[1].size()+1);
-
+
cells.resize (step_sizes[2].size()*step_sizes[1].size()*step_sizes[0].size());
for (unsigned int z=0; z<step_sizes[2].size(); ++z)
for (unsigned int y=0; y<step_sizes[1].size(); ++y)
cells[c].material_id = 0;
}
break;
-
+
}
default:
Assert (false, ExcNotImplemented());
}
- tria.create_triangulation (points, cells, SubCellData());
+ tria.create_triangulation (points, cells, SubCellData());
if (colorize)
{
// faces of all cells and set
// boundary indicator to the
// correct value if it was 0.
-
+
// use a large epsilon to
// compare numbers to avoid
// roundoff problems.
*std::min_element (step_sizes[i].begin(),
step_sizes[i].end()));
const double epsilon = 0.01 * min_size;
-
+
// actual code is external since
// 1-D is different from 2/3D.
colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon);
const bool colorize)
{
// contributed by Yaqi Wang 2006
- Assert(spacing.size() == 1,
+ Assert(spacing.size() == 1,
ExcInvalidRepetitionsDimension(1));
const unsigned int n_cells = material_id.size(0);
- Assert(spacing[0].size() == n_cells,
+ Assert(spacing[0].size() == n_cells,
ExcInvalidRepetitionsDimension(1));
double delta = std::numeric_limits<double>::max();
Assert (spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
delta = std::min (delta, spacing[0][i]);
}
-
+
// generate the necessary points
std::vector<Point<1> > points;
double ax = p[0];
SubCellData t;
GridTools::delete_unused_vertices (points, cells, t);
- tria.create_triangulation (points, cells, t);
+ tria.create_triangulation (points, cells, t);
// set boundary indicator
if (colorize)
const bool colorize)
{
// contributed by Yaqi Wang 2006
- Assert(spacing.size() == 2,
+ Assert(spacing.size() == 2,
ExcInvalidRepetitionsDimension(2));
std::vector<unsigned int> repetitions(2);
Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
delta = std::min (delta, spacing[i][j]);
}
- Assert(material_id.size(i) == repetitions[i],
+ Assert(material_id.size(i) == repetitions[i],
ExcInvalidRepetitionsDimension(i));
}
-
+
// generate the necessary points
std::vector<Point<2> > points;
double ay = p[1];
SubCellData t;
GridTools::delete_unused_vertices (points, cells, t);
- tria.create_triangulation (points, cells, t);
+ tria.create_triangulation (points, cells, t);
// set boundary indicator
if (colorize)
{
// contributed by Yaqi Wang 2006
const unsigned int dim = 3;
-
- Assert(spacing.size() == dim,
+
+ Assert(spacing.size() == dim,
ExcInvalidRepetitionsDimension(dim));
std::vector<unsigned int> repetitions(dim);
Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
delta = std::min (delta, spacing[i][j]);
}
- Assert(material_id.size(i) == repetitions[i],
+ Assert(material_id.size(i) == repetitions[i],
ExcInvalidRepetitionsDimension(i));
}
SubCellData t;
GridTools::delete_unused_vertices (points, cells, t);
- tria.create_triangulation (points, cells, t);
+ tria.create_triangulation (points, cells, t);
// set boundary indicator
if (colorize)
const Point<dim> &p2,
const double epsilon)
{
-
+
// run through all faces and check
// if one of their center coordinates matches
// one of the corner points. Comparisons
// are made using an epsilon which
// should be smaller than the smallest cell
// diameter.
-
+
typename Triangulation<dim>::raw_face_iterator face = tria.begin_raw_face(),
endface = tria.end_face();
for (; face!=endface; ++face)
// but we could not find
// on which boundary.
Assert (false, ExcInternalError());
-
+
}
}
for (typename Triangulation<dim>::cell_iterator cell = tria.begin();
+// Implementation for 1D only
+template <int dim>
+void GridGenerator::truncated_cone (Triangulation<dim> &,
+ const double,
+ const double,
+ const double)
+{
+ Assert (false, ExcNotImplemented());
+}
+
+
+
// Implementation for 1D only
template <int dim>
void GridGenerator::hyper_shell (Triangulation<dim> &,
1, 0, 2,
9, 8,10
};
-
+
std::vector<CellData<dim> > cells(9);
k = 0;
for (unsigned int i0=0;i0<3;++i0)
}
+
+template <int dim>
+void GridGenerator::truncated_cone (Triangulation<dim> &triangulation,
+ const double radius_0,
+ const double radius_1,
+ const double half_length)
+{
+ Point<dim> vertices_tmp[4];
+
+ vertices_tmp[0] = Point<dim> (-half_length, -radius_0);
+ vertices_tmp[1] = Point<dim> (half_length, -radius_1);
+ vertices_tmp[2] = Point<dim> (-half_length, radius_0);
+ vertices_tmp[3] = Point<dim> (half_length, radius_1);
+
+ const std::vector<Point<dim> > vertices (&vertices_tmp[0], &vertices_tmp[4]);
+ unsigned int cell_vertices[1][GeometryInfo<dim>::vertices_per_cell];
+
+ for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
+ cell_vertices[0][i] = i;
+
+ std::vector<CellData<dim> > cells (1, CellData<dim> ());
+
+ for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
+ cells[0].vertices[i] = cell_vertices[0][i];
+
+ cells[0].material_id = 0;
+ triangulation.create_triangulation (vertices, cells, SubCellData ());
+
+ typename Triangulation<dim>::cell_iterator cell = triangulation.begin ();
+
+ cell->face (0)->set_boundary_indicator (1);
+ cell->face (1)->set_boundary_indicator (2);
+
+ for (unsigned int i = 2; i < 4; ++i)
+ cell->face (i)->set_boundary_indicator (0);
+}
+
+
+
//TODO: Colorize edges as circumference, left and right radius
// Implementation for 2D only
template <int dim>
{3, 4, 6, 7}};
std::vector<CellData<dim> > cells (3, CellData<dim>());
-
- for (unsigned int i=0; i<3; ++i)
+
+ for (unsigned int i=0; i<3; ++i)
{
for (unsigned int j=0; j<4; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
};
-
+
tria.create_triangulation (
std::vector<Point<dim> >(&vertices[0], &vertices[8]),
cells,
p+Point<dim>(+1,+1)*(radius/std::sqrt(2.0)*a),
p+Point<dim>(-1,+1)*(radius/std::sqrt(2.0)),
p+Point<dim>(+1,+1)*(radius/std::sqrt(2.0)) };
-
+
const int cell_vertices[5][4] = {{0, 1, 2, 3},
{0, 2, 6, 4},
{2, 3, 4, 5},
{6, 4, 7, 5}};
std::vector<CellData<dim> > cells (5, CellData<dim>());
-
- for (unsigned int i=0; i<5; ++i)
+
+ for (unsigned int i=0; i<5; ++i)
{
for (unsigned int j=0; j<4; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
};
-
+
tria.create_triangulation (
std::vector<Point<dim> >(&vertices[0], &vertices[8]),
cells,
{
Assert ((inner_radius > 0) && (inner_radius < outer_radius),
ExcInvalidRadii ());
-
+
const double pi = numbers::PI;
-
+
// determine the number of cells
// for the grid. if not provided by
// the user determine it such that
};
std::vector<CellData<dim> > cells (N, CellData<dim>());
-
- for (unsigned int i=0; i<N; ++i)
+
+ for (unsigned int i=0; i<N; ++i)
{
cells[i].vertices[0] = i;
cells[i].vertices[1] = (i+1)%N;
cells[i].vertices[2] = N+i;
cells[i].vertices[3] = N+((i+1)%N);
-
+
cells[i].material_id = 0;
};
-
+
tria.create_triangulation (
vertices, cells, SubCellData());
break;
default:
f->set_boundary_indicator(0);
- break;
+ break;
}
++f;
}
p+Point<dim>(+1,+1)*(radius/std::sqrt(2.0)*a),
p+Point<dim>(0,+1)*radius,
p+Point<dim>(+1,+1)*(radius/std::sqrt(2.0)) };
-
+
const int cell_vertices[5][4] = {{0, 1, 2, 3},
{2, 3, 4, 5},
{1, 7, 3, 5},
{6, 4, 7, 5}};
std::vector<CellData<dim> > cells (4, CellData<dim>());
-
- for (unsigned int i=0; i<4; ++i)
+
+ for (unsigned int i=0; i<4; ++i)
{
for (unsigned int j=0; j<4; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
};
-
+
tria.create_triangulation (
std::vector<Point<dim> >(&vertices[0], &vertices[8]),
cells,
{
Assert ((inner_radius > 0) && (inner_radius < outer_radius),
ExcInvalidRadii ());
-
+
const double pi = numbers::PI;
// determine the number of cells
// for the grid. if not provided by
std::vector<CellData<dim> > cells (N, CellData<dim>());
-
- for (unsigned int i=0; i<N; ++i)
+
+ for (unsigned int i=0; i<N; ++i)
{
cells[i].vertices[0] = i;
cells[i].vertices[1] = (i+1)%(N+1);
cells[i].vertices[2] = N+1+i;
cells[i].vertices[3] = N+1+((i+1)%(N+1));
-
+
cells[i].material_id = 0;
};
-
+
tria.create_triangulation (vertices, cells, SubCellData());
}
{
const double rl2=(right+left)/2;
const double len = (right-left)/2.;
-
+
const Point<dim> vertices[20] = {
Point<dim>(left, left , -len/2.),
Point<dim>(rl2, left , -len/2.),
33,32,34,
41,40,42
};
-
+
std::vector<CellData<dim> > cells(27);
k = 0;
for (unsigned int z=0;z<3;++z)
+template <int dim>
+void GridGenerator::truncated_cone (Triangulation<dim> &triangulation,
+ const double radius_0,
+ const double radius_1,
+ const double half_length)
+{
+ Assert (dim <= 3, ExcNotImplemented());
+ // Determine number of cells and vertices
+ const unsigned int n_cells = std::floor (half_length / std::max (radius_0, radius_1) + 0.5);
+ const unsigned int n_vertices = 4 * (n_cells + 1);
+ std::vector<Point<dim> > vertices_tmp(n_vertices);
+
+ vertices_tmp[0] = Point<dim> (-half_length, 0, -radius_0);
+ vertices_tmp[1] = Point<dim> (-half_length, radius_0, 0);
+ vertices_tmp[2] = Point<dim> (-half_length, -radius_0, 0);
+ vertices_tmp[3] = Point<dim> (-half_length, 0, radius_0);
+
+ const double dx = 2 * half_length / n_cells;
+
+ for (unsigned int i = 0; i < n_cells; ++i) {
+ vertices_tmp[4 * (i + 1)]
+ = vertices_tmp[4 * i] +
+ Point<dim> (dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length);
+ vertices_tmp[4 * i + 5]
+ = vertices_tmp[4 * i + 1] +
+ Point<dim> (dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0);
+ vertices_tmp[4 * i + 6]
+ = vertices_tmp[4 * i + 2] +
+ Point<dim> (dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0);
+ vertices_tmp[4 * i + 7]
+ = vertices_tmp[4 * i + 3] +
+ Point<dim> (dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length);
+ }
+
+ const std::vector<Point<dim> > vertices (&vertices_tmp[0], &vertices_tmp[n_vertices]);
+ Table<2,unsigned int> cell_vertices(n_cells,GeometryInfo<dim>::vertices_per_cell);
+
+ for (unsigned int i = 0; i < n_cells; ++i)
+ for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j)
+ cell_vertices[i][j] = 4 * i + j;
+
+ std::vector<CellData<dim> > cells (n_cells, CellData<dim> ());
+
+ for (unsigned int i = 0; i < n_cells; ++i) {
+ for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j)
+ cells[i].vertices[j] = cell_vertices[i][j];
+
+ cells[i].material_id = 0;
+ }
+
+ triangulation.create_triangulation (vertices, cells, SubCellData ());
+
+ for (typename Triangulation<dim>::cell_iterator cell = triangulation.begin ();
+ cell != triangulation.end (); ++cell) {
+ if (cell->vertex (0) (0) == -half_length) {
+ cell->face (4)->set_boundary_indicator (1);
+
+ for (unsigned int i = 0; i < 4; ++i)
+ cell->line (i)->set_boundary_indicator (0);
+ }
+
+ if (cell->vertex (4) (0) == half_length) {
+ cell->face (5)->set_boundary_indicator (2);
+
+ for (unsigned int i = 4; i < 8; ++i)
+ cell->line (i)->set_boundary_indicator (0);
+ }
+
+ for (unsigned int i = 0; i < 4; ++i)
+ cell->face (i)->set_boundary_indicator (0);
+ }
+}
+
+
// Implementation for 3D only
template <int dim>
void
{12, 13, 21, 22, 15, 16, 24, 25}};
std::vector<CellData<dim> > cells (7, CellData<dim>());
-
- for (unsigned int i=0; i<7; ++i)
+
+ for (unsigned int i=0; i<7; ++i)
{
for (unsigned int j=0; j<8; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
{8, 0, 12, 4, 11, 3, 15, 7}, // left
{8, 9, 0, 1, 11, 10, 3, 2}, // front
{12, 4, 13, 5, 15, 7, 14, 6}}; // back
-
+
std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
-
- for (unsigned int i=0; i<n_cells; ++i)
+
+ for (unsigned int i=0; i<n_cells; ++i)
{
for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
const double half_length)
{
Assert (dim <= 3, ExcNotImplemented());
-
+
// Copy the base from hyper_ball<dim>
// and transform it to yz
const double d = radius/std::sqrt(2.0);
vertices[i](1) = -vertices[i](0);
vertices[i](0) = h;
}
-
+
int cell_vertices[10][8] = {
{0, 1, 8, 9, 2, 3, 10, 11},
{0, 2, 8, 10, 6, 4, 14, 12},
for (unsigned int i=0;i<5;++i)
for (unsigned int j=0;j<8;++j)
cell_vertices[i+5][j] = cell_vertices[i][j]+8;
-
+
std::vector<CellData<dim> > cells (10, CellData<dim>());
-
- for (unsigned int i=0; i<10; ++i)
+
+ for (unsigned int i=0; i<10; ++i)
{
for (unsigned int j=0; j<8; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
// above
typename Triangulation<dim>::cell_iterator cell = tria.begin();
typename Triangulation<dim>::cell_iterator end = tria.end();
-
+
for (; cell != end; ++cell)
for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
if (cell->at_boundary(i))
std::vector<CellData<dim> > cells (6, CellData<dim>());
- for (unsigned int i=0; i<6; ++i)
+ for (unsigned int i=0; i<6; ++i)
{
for (unsigned int j=0; j<8; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
// element
if (cell->face(i)->center()(0) < center(0)+1.e-5) {
cell->face(i)->set_boundary_indicator(1);
- for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j)
+ for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j)
cell->face(i)->line(j)->set_boundary_indicator(1);
}
}
// If the center is not on the plane x=0, this is a curvilinear
// element
if (cell->face(i)->center()(0) > center(0)+1.e-5) {
- for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j)
+ for (unsigned int j=0;j<GeometryInfo<dim>::lines_per_face;++j)
cell->face(i)->line(j)->set_boundary_indicator(0);
}
}
const double orad = outer_radius/std::sqrt(3.0);
std::vector<Point<dim> > vertices;
std::vector<CellData<dim> > cells;
-
+
// Start with the shell bounded by
// two nested cubes
if (n <= 6)
vertices.push_back(p+hexahedron[i]*irad);
for (unsigned int i=0;i<8;++i)
vertices.push_back(p+hexahedron[i]*orad);
-
+
const unsigned int n_cells = 6;
const int cell_vertices[n_cells][8] =
{{8, 9, 10, 11, 0, 1, 2, 3}, // bottom
{8, 0, 10, 2, 12, 4, 14, 6}, // left
{8, 9, 0, 1, 12, 13, 4, 5}, // front
{10, 2, 11, 3, 14, 6, 15, 7}}; // back
-
+
cells.resize(n_cells, CellData<dim>());
-
- for (unsigned int i=0; i<n_cells; ++i)
+
+ for (unsigned int i=0; i<n_cells; ++i)
{
for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
{ 2, 8, 11, 6},
{ 12, 0, 2, 8},
{ 1, 10, 12, 0}};
-
+
cells.resize(n_cells, CellData<dim>());
-
- for (unsigned int i=0; i<n_cells; ++i)
+
+ for (unsigned int i=0; i<n_cells; ++i)
{
for (unsigned int j=0; j<4; ++j)
{
{
Assert(false, ExcIndexRange(n, 1, 7));
}
-
+
tria.create_triangulation (vertices, cells,
SubCellData()); // no boundary
// information
-
+
if (colorize)
colorize_hyper_shell(tria, p, inner_radius, outer_radius);
}
(++cell)->face(2)->set_boundary_indicator(1);
(++cell)->face(0)->set_boundary_indicator(1);
(++cell)->face(2)->set_boundary_indicator(1);
- (++cell)->face(0)->set_boundary_indicator(1);
+ (++cell)->face(0)->set_boundary_indicator(1);
}
else
// For higher polyhedra, this is regular.
std::vector<CellData<dim> > cells (5, CellData<dim>());
- for (unsigned int i=0; i<5; ++i)
+ for (unsigned int i=0; i<5; ++i)
{
for (unsigned int j=0; j<8; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
{
Assert ((inner_radius > 0) && (inner_radius < outer_radius),
ExcInvalidRadii ());
-
+
const double pi = numbers::PI;
-
+
// determine the number of cells
// for the grid. if not provided by
// the user determine it such that
j*length/N_z);
vertices_3d.push_back (v);
}
-
+
std::vector<CellData<dim> > cells (N_r*N_z, CellData<dim>());
-
+
for (unsigned int j=0; j<N_z; ++j)
- for (unsigned int i=0; i<N_r; ++i)
+ for (unsigned int i=0; i<N_r; ++i)
{
cells[i+j*N_r].vertices[0] = i + (j+1)*2*N_r;
cells[i+j*N_r].vertices[1] = (i+1)%N_r + (j+1)*2*N_r;
cells[i+j*N_r].vertices[5] = N_r+((i+1)%N_r) + (j+1)*2*N_r;
cells[i+j*N_r].vertices[6] = N_r+i + j*2*N_r;
cells[i+j*N_r].vertices[7] = N_r+((i+1)%N_r) + j*2*N_r;
-
+
cells[i+j*N_r].material_id = 0;
}
-
+
tria.create_triangulation (
vertices_3d, cells, SubCellData());
}
PreconditionJacobi<SparseMatrix<double> > prec;
prec.initialize(S, 1.2);
FilteredMatrix<Vector<double> > PF (prec);
-
+
SolverControl control (n_dofs, 1.e-10, false, false);
GrowingVectorMemory<Vector<double> > mem;
SolverCG<Vector<double> > solver (control, mem);
-
+
Vector<double> f(n_dofs);
-
+
SF.add_constraints(m);
SF.apply_constraints (f, true);
solver.solve(SF, u, f, PF);
{
// first provide everything that is
// needed for solving a Laplace
- // equation.
- MappingQ1<dim> mapping_q1;
+ // equation.
+ MappingQ1<dim> mapping_q1;
FE_Q<dim> q1(1);
-
+
DoFHandler<dim> dof_handler(tria);
dof_handler.distribute_dofs(q1);
SparsityPattern sparsity_pattern (dof_handler.n_dofs (), dof_handler.n_dofs (),
dof_handler.max_couplings_between_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
sparsity_pattern.compress ();
-
+
SparseMatrix<double> S(sparsity_pattern);
-
+
QGauss4<dim> quadrature;
-
+
MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
// set up the boundary values for
}
}
}
-
+
// solve the dim problems with
// different right hand sides.
Vector<double> us[dim];
for (unsigned int i=0; i<dim; ++i)
us[i].reinit (dof_handler.n_dofs());
-
+
// solve linear systems in parallel
Threads::TaskGroup<> tasks;
for (unsigned int i=0; i<dim; ++i)
tasks += Threads::new_task (&GridGenerator::laplace_solve,
S, m[i], us[i]);
tasks.join_all ();
-
+
// change the coordinates of the
// points of the triangulation
// according to the computed values
#if deal_II_dimension == 2
template<int dim>
-void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
+void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
const double inner_radius,
const double outer_radius,
const double, // width,
}
}
}
- }
+ }
}
double eps = 1e-3 * outer_radius;
cell = triangulation.begin_active();
cell->face(f)->set_boundary_indicator(2);
else if(std::abs(dy - outer_radius) < eps)
cell->face(f)->set_boundary_indicator(3);
- else
+ else
cell->face(f)->set_boundary_indicator(4);
} else {
double d = (cell->face(f)->center() - center).norm();
cell->face(f)->set_boundary_indicator(0);
}
}
- }
+ }
}
#endif
#if deal_II_dimension == 3
template<int dim>
-void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangulation,
+void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangulation,
const double inner_radius,
const double outer_radius,
const double L,
{
Assert(inner_radius < outer_radius,
ExcMessage("outer_radius has to be bigger than inner_radius."));
- Assert(L > 0,
+ Assert(L > 0,
ExcMessage("Must give positive extension L"));
Assert(Nz >= 1, ExcLowerRange(1, Nz));
-
+
GridGenerator::cylinder_shell (triangulation,
L, inner_radius, outer_radius,
8,
Nz);
-
+
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
}
}
}
- }
+ }
}
double eps = 1e-3 * outer_radius;
cell = triangulation.begin_active();
double dx = cell->face(f)->center()(0);
double dy = cell->face(f)->center()(1);
double dz = cell->face(f)->center()(2);
-
+
if(colorize) {
if(std::abs(dx + outer_radius) < eps)
cell->face(f)->set_boundary_indicator(0);
-
+
else if(std::abs(dx - outer_radius) < eps)
cell->face(f)->set_boundary_indicator(1);
-
+
else if(std::abs(dy + outer_radius) < eps)
cell->face(f)->set_boundary_indicator(2);
-
+
else if(std::abs(dy - outer_radius) < eps)
cell->face(f)->set_boundary_indicator(3);
-
+
else if(std::abs(dz) < eps)
cell->face(f)->set_boundary_indicator(4);
else if(std::abs(dz - L) < eps)
cell->face(f)->set_boundary_indicator(5);
-
+
else {
cell->face(f)->set_boundary_indicator(6);
for(unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
cell->face(f)->set_boundary_indicator(0);
}
}
- }
+ }
}
#endif
const bool);
template void
-GridGenerator::subdivided_hyper_rectangle<deal_II_dimension>
+GridGenerator::subdivided_hyper_rectangle<deal_II_dimension>
(Triangulation<deal_II_dimension> &,
const std::vector<unsigned int>&,
const Point<deal_II_dimension>&,
const Point<deal_II_dimension>&, bool);
template void
-GridGenerator::subdivided_hyper_rectangle<deal_II_dimension>
+GridGenerator::subdivided_hyper_rectangle<deal_II_dimension>
(Triangulation<deal_II_dimension> &,
const std::vector<std::vector<double> >&,
const Point<deal_II_dimension>&,
Triangulation<deal_II_dimension>&, double, double);
+template void
+GridGenerator::truncated_cone (
+ Triangulation<deal_II_dimension>&, double, double, double);
+
template void
GridGenerator::hyper_L (
Triangulation<deal_II_dimension>&, double, double);
const Point<deal_II_dimension>&, double, double, unsigned int);
-template void
+template void
GridGenerator::hyper_cube_with_cylindrical_hole (
Triangulation<deal_II_dimension> &,
const double, const double, const double, const unsigned int, bool);
get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
{
const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
-
+
// same algorithm as above
const unsigned int dim = 3;
const Point<dim> vector_from_axis = (middle-point_on_axis) -
get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
}
-
+
template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_between_points (
{
unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
Assert(points.size()==m*m, ExcInternalError());
-
+
std::vector<Point<3> > lp0(m);
std::vector<Point<3> > lp1(m);
get_intermediate_points_on_line(quad->line(0), lp0);
get_intermediate_points_on_line(quad->line(1), lp1);
-
+
std::vector<Point<3> > lps(m);
for (unsigned int i=0; i<m; ++i)
{
get_intermediate_points_between_points(lp0[i], lp1[i], lps);
-
+
for (unsigned int j=0; j<m; ++j)
points[i*m+j]=lps[j];
}
const Point<dim> vector_from_axis = (vertex-point_on_axis) -
((vertex-point_on_axis) * direction) * direction;
-
+
face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm());
}
}
template <int dim>
double
-CylinderBoundary<dim>::get_radius () const
+CylinderBoundary<dim>::get_radius () const
{
return radius;
}
+//======================================================================//
+
+template<int dim>
+ConeBoundary<dim>::ConeBoundary (const double radius_0,
+ const double radius_1,
+ const Point<dim> x_0,
+ const Point<dim> x_1)
+ :
+ radius_0 (radius_0),
+ radius_1 (radius_1),
+ x_0 (x_0),
+ x_1 (x_1)
+{}
+
+
+
+template<int dim>
+double ConeBoundary<dim>::get_radius (Point<dim> x) const
+{
+ for (unsigned int i = 0; i < dim; ++i)
+ if ((x_1 (i) - x_0 (i)) != 0)
+ return (radius_1 - radius_0) * (x (i) - x_0 (i)) / (x_1 (i) - x_0 (i)) + radius_0;
+
+ return 0;
+}
+
+
+
+template<int dim>
+void
+ConeBoundary<dim>::
+get_intermediate_points_between_points (const Point<dim> &p0,
+ const Point<dim> &p1,
+ std::vector<Point<dim> > &points) const
+{
+ const unsigned int n = points.size ();
+ const Point<dim> axis = x_1 - x_0;
+ const Point<dim> dx = (p1 - p0) / (n + 1);
+
+ Assert (n > 0, ExcInternalError ());
+
+ for (unsigned int i = 0; i < n; ++i)
+ {
+ // Compute the current point.
+ const Point<dim> x_i = p0 + (i + 1) * dx;
+ // To project this point on the
+ // boundary of the cone we first
+ // compute the orthogonal
+ // projection of this point onto
+ // the axis of the cone.
+ const double c = (x_i - x_0) * axis / axis.square ();
+ const Point<dim> x_ip = x_0 + c * axis;
+ // If the middle point is located
+ // on the axis hand over the
+ // current point itself, else
+ // hand over the projection of it
+ // on the boundary.
+ if (x_ip.norm () <= 1e-10 * x_i.norm ())
+ points[i] = x_i;
+ else {
+ // Compute the projection of
+ // the middle point on the
+ // boundary of the cone.
+ points[i] = x_ip + get_radius (x_ip) * (x_i - x_ip) / (x_i - x_ip).norm ();
+ }
+ }
+}
+
+template<int dim>
+Point<dim>
+ConeBoundary<dim>::
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+ const Point<dim> axis = x_1 - x_0;
+ // Compute the middle point of the line.
+ const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+ // To project it on the boundary of
+ // the cone we first compute the
+ // orthogonal projection of the
+ // middle point onto the axis of
+ // the cone.
+ const double c = (middle - x_0) * axis / axis.square ();
+ const Point<dim> middle_p = x_0 + c * axis;
+ // If the middle point is located
+ // on the axis return the middle
+ // point itself, else return the
+ // projection of it on the
+ // boundary.
+ if (middle_p.norm () <= 1e-10 * middle.norm ())
+ return middle;
+ else {
+ // Compute the projection of the
+ // middle point on the boundary
+ // of the cone.
+ return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
+ }
+}
+
+
+
+template<>
+Point<3>
+ConeBoundary<3>::get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
+{
+ const Point<3> axis = x_1 - x_0;
+ // Compute the middle point of the
+ // quad.
+ const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+ // Same algorithm as above: To
+ // project it on the boundary of
+ // the cone we first compute the
+ // orthogonal projection of the
+ // middle point onto the axis of
+ // the cone.
+ const double c = (middle - x_0) * axis / axis.square ();
+ const Point<3> middle_p = x_0 + c * axis;
+ // If the middle point is located
+ // on the axis return the middle
+ // point itself, else return the
+ // projection of it on the
+ // boundary.
+ if (middle_p.norm () <= 1e-10 * middle.norm ())
+ return middle;
+ else {
+ // Compute the projection of the
+ // middle point on the boundary
+ // of the cone.
+ return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
+ }
+}
+
+
+
+template<int dim>
+Point<dim>
+ConeBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
+{
+ Assert (false, ExcImpossibleInDim (dim));
+
+ return Point<dim>();
+}
+
+
+
+template<int dim>
+void
+ConeBoundary<dim>::
+get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+ std::vector<Point<dim> > &points) const
+{
+ if (points.size () == 1)
+ points[0] = get_new_point_on_line (line);
+ else
+ get_intermediate_points_between_points (line->vertex (0), line->vertex (1), points);
+}
+
+
+#if deal_II_dimension == 3
+
+template<>
+void
+ConeBoundary<3>::
+get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
+ std::vector<Point<3> > &points) const
+{
+ if (points.size () == 1)
+ points[0] = get_new_point_on_quad (quad);
+ else
+ {
+ unsigned int n = static_cast<unsigned int> (std::sqrt (static_cast<double> (points.size ())));
+
+ Assert (points.size () == n * n, ExcInternalError ());
+
+ std::vector<Point<3> > points_line_0 (n);
+ std::vector<Point<3> > points_line_1 (n);
+
+ get_intermediate_points_on_line (quad->line (0), points_line_0);
+ get_intermediate_points_on_line (quad->line (1), points_line_1);
+
+ std::vector<Point<3> > points_line_segment (n);
+
+ for (unsigned int i = 0; i < n; ++i)
+ {
+ get_intermediate_points_between_points (points_line_0[i],
+ points_line_1[i],
+ points_line_segment);
+
+ for (unsigned int j = 0; j < n; ++j)
+ points[i * n + j] = points_line_segment[j];
+ }
+ }
+}
+
+#endif
+
+template <int dim>
+void
+ConeBoundary<dim>::
+get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
+ std::vector<Point<dim> > &) const
+{
+ Assert (false, ExcImpossibleInDim (dim));
+}
+
+
+#if deal_II_dimension == 1
+
+template<>
+void
+ConeBoundary<1>::
+get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+ Boundary<1>::FaceVertexNormals &) const
+{
+ Assert (false, ExcImpossibleInDim (1));
+}
+
+#endif
+
+template<int dim>
+void
+ConeBoundary<dim>::
+get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+{
+ const Point<dim> axis = x_1 - x_0;
+
+ for (unsigned int vertex = 0; vertex < GeometryInfo<dim>::vertices_per_cell; ++vertex)
+ {
+ // Compute the orthogonal
+ // projection of the vertex onto
+ // the axis of the cone.
+ const double c = (face->vertex (vertex) - x_0) * axis / axis.square ();
+ const Point<dim> vertex_p = x_0 + c * axis;
+ // Then compute the vector
+ // pointing from the point
+ // <tt>vertex_p</tt> on the axis
+ // to the vertex.
+ const Point<dim> axis_to_vertex = face->vertex (vertex) - vertex_p;
+
+ face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm ();
+ }
+}
+
+
//======================================================================//
template <int dim, int spacedim>
HyperBallBoundary<dim,spacedim>::get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
{
Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_line (line);
-
+
middle -= center;
double r=0;
r=radius;
// project to boundary
middle *= r / std::sqrt(middle.square());
-
+
middle += center;
return middle;
}
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const
{
Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_quad (quad);
-
+
middle -= center;
double r=0;
r=radius;
// project to boundary
middle *= r / std::sqrt(middle.square());
-
+
middle += center;
return middle;
}
get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
}
-
+
template <int dim, int spacedim>
void
HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
const Point<spacedim> v0=p0-center,
v1=p1-center;
const double length=std::sqrt((v1-v0).square());
-
+
double eps=1e-12;
double r=0;
if (compute_radius_automatically)
else
r=radius;
-
+
const double r2=r*r;
Assert(std::fabs(v0.square()-r2)<eps*r2, ExcInternalError());
Assert(std::fabs(v1.square()-r2)<eps*r2, ExcInternalError());
-
+
const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square()));
const double d_alpha=alpha/(n+1);
const Point<spacedim> pm=0.5*(v0+v1);
-
+
const double h=std::sqrt(pm.square());
-
+
double beta=0;
unsigned int left_index=0, right_index=0;
if ((n+1)%2==0)
right_index=n/2;
left_index=n/2-1;
}
-
+
// n even: m=n/2,
// n odd: m=(n-1)/2
const unsigned int m=n/2;
{
Assert(right_index<n, ExcInternalError());
Assert(left_index+1>=1, ExcInternalError());
-
+
double d=h*std::tan(beta);
points[right_index]=pm+d/length*(v1-v0);
points[left_index]=pm-d/length*(v1-v0);
}
-
+
// project the points from the
// straight line to the
{
unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
Assert(points.size()==m*m, ExcInternalError());
-
+
std::vector<Point<3> > lp0(m);
std::vector<Point<3> > lp1(m);
-
+
get_intermediate_points_on_line(quad->line(0), lp0);
get_intermediate_points_on_line(quad->line(1), lp1);
-
+
std::vector<Point<3> > lps(m);
for (unsigned int i=0; i<m; ++i)
{
get_intermediate_points_between_points(lp0[i], lp1[i], lps);
-
+
for (unsigned int j=0; j<m; ++j)
points[i*m+j]=lps[j];
}
template <int dim, int spacedim>
Point<spacedim>
-HyperBallBoundary<dim,spacedim>::get_center () const
+HyperBallBoundary<dim,spacedim>::get_center () const
{
return center;
}
template <int dim, int spacedim>
double
-HyperBallBoundary<dim,spacedim>::get_radius () const
+HyperBallBoundary<dim,spacedim>::get_radius () const
{
Assert(!compute_radius_automatically, ExcRadiusNotSet());
return radius;
template <int dim>
Point<dim>
HalfHyperBallBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
{
// check whether center of object is
// at x==0, since then it belongs
template <>
Point<1>
HalfHyperBallBoundary<1>::
-get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
+get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
{
Assert (false, ExcInternalError());
return Point<1>();
template <int dim>
Point<dim>
HalfHyperBallBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
{
const Point<dim> quad_center = quad->center();
if (quad_center(0) == this->center(0))
template <int dim>
Point<dim>
HalfHyperShellBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
{
switch (dim)
{
// as in the base class
return HyperShellBoundary<dim>::get_new_point_on_line (line);
}
-
+
// in 3d, a line is a straight
// line if it is on the symmetry
// plane and if not both of its
// inner or outer sphere
case 3:
{
-
+
if (((line->vertex(0)(0) == this->center(0))
&&
(line->vertex(1)(0) == this->center(0)))
default:
Assert (false, ExcNotImplemented());
}
-
+
return Point<dim>();
}
quad->vertex(2) + quad->vertex(3) )/4;
const Point<dim> quad_center_offset = quad_center - this->center;
-
+
if (std::fabs (quad->line(0)->center().distance(this->center) -
quad->line(1)->center().distance(this->center))
< 1e-12 * outer_radius)
// lines 0 and 1 are radial
const double needed_radius
= quad->line(0)->center().distance(this->center);
-
+
return (this->center +
quad_center_offset/quad_center_offset.norm() * needed_radius);
}
// lines 2 and 3 are radial
const double needed_radius
= quad->line(2)->center().distance(this->center);
-
+
return (this->center +
quad_center_offset/quad_center_offset.norm() * needed_radius);
}
else
Assert (false, ExcInternalError());
}
-
+
// otherwise we are on the outer or
// inner part of the shell. proceed
// as in the base class
// as in the base class
HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
}
-
+
// in 3d, a line is a straight
// line if it is on the symmetry
// plane and if not both of its
// inner or outer sphere
case 3:
{
-
+
if (((line->vertex(0)(0) == this->center(0))
&&
(line->vertex(1)(0) == this->center(0)))
std::vector<Point<dim> > &points) const
{
Assert (dim < 3, ExcNotImplemented());
-
+
// check whether center of object is
// at x==0, since then it belongs
// to the plane part of the
// explicit instantiations
-template class HyperBallBoundary<deal_II_dimension>;
template class CylinderBoundary<deal_II_dimension>;
+template class ConeBoundary<deal_II_dimension>;
+template class HyperBallBoundary<deal_II_dimension>;
template class HalfHyperBallBoundary<deal_II_dimension>;
template class HyperShellBoundary<deal_II_dimension>;
template class HalfHyperShellBoundary<deal_II_dimension>;