const double b = a/2.0;
const double c = d/2.0;
// And so are these
- const double hb = a*std::sqrt(3.0)/2.0;
- const double hc = d*std::sqrt(3.0)/2.0;
+ const double hb = radius*std::sqrt(3.0)/4.0;
+ const double hc = radius*std::sqrt(3.0)/2.0;
Point<dim> vertices[16] = {
center+Point<dim>( 0, d, -d),
if (cell->face(i)->boundary_indicator() == 255)
continue;
- // If the center is on the plane x=0, this is a plane
+ // If the center is on the plane x=0, this is a planar
// element
- if (cell->face(i)->center()(0) < center(0)+1.e-5)
+ 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)
+ cell->face(i)->line(j)->set_boundary_indicator(1);
+ }
+ }
+ // With this loop we restore back the indicator of the outer lines
+ for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+ {
+ if (cell->face(i)->boundary_indicator() == 255)
+ continue;
+
+ // 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)
+ cell->face(i)->line(j)->set_boundary_indicator(0);
+ }
}
++cell;
}
}
-
// Implementation for 3D only
template <int dim>
void GridGenerator::hyper_shell (Triangulation<dim>&,
2.00000 -1.87868 0.878680
2.00000 1.12132 2.12132
2.00000 -3.12132 2.12132
- 3.83712 0.0606602 -1.06066
- 3.83712 -2.06066 -1.06066
- 2.76096 -0.560660 -0.439340
- 2.76096 -1.43934 -0.439340
- 2.76096 -0.560660 0.439340
- 2.76096 -1.43934 0.439340
- 3.83712 0.0606602 1.06066
- 3.83712 -2.06066 1.06066
+ 4.59808 0.0606602 -1.06066
+ 4.59808 -2.06066 -1.06066
+ 3.29904 -0.560660 -0.439340
+ 3.29904 -1.43934 -0.439340
+ 3.29904 -0.560660 0.439340
+ 3.29904 -1.43934 0.439340
+ 4.59808 0.0606602 1.06066
+ 4.59808 -2.06066 1.06066
object "cells" class array type int rank 1 shape 8 items 6 data follows
0 2 8 10 1 3 9 11
0 6 8 14 2 4 10 12