#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary_lib.h>
+#include <grid/grid_reordering.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <hp/fe_values.h>
RightHandSide<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
- double product = 1;
- for (unsigned int d=0; d<dim; ++d)
- product *= (p[d]+1);
- return product;
+ switch (dim)
+ {
+ case 2:
+ {
+ double product = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ product *= (p[d]+1);
+ return product;
+ }
+
+ case 3:
+ return (p[0]>std::fabs(p[1]) ? 1 : 0);
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ return 0.;
}
template <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
dof_handler (triangulation),
- max_degree (dim == 2 ? 7 : 4)
+ max_degree (dim == 2 ? 7 : 5)
{
for (unsigned int degree=2; degree<=max_degree; ++degree)
{
const std::string filename = "solution-" +
Utilities::int_to_string (cycle, 2) +
- ".vtk";
+ ".gmv";
DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.build_patches ();
std::ofstream output (filename.c_str());
- data_out.write_vtk (output);
+ data_out.write_gmv (output);
}
}
-template <>
-void LaplaceProblem<3>::create_coarse_grid ()
+namespace BreastPhantom
{
- const unsigned int dim = 3;
-
- // GridGenerator::hyper_cube (triangulation);
- // triangulation.refine_global (1);
-
- // Create a hollow cube, in analogy to the 2D example.
- // The grid generation is done in two steps. First the
- // cell data is created on a uniform grid. In the
- // second step, the unused vertices are removed.
- const unsigned char hollow [4][4] = {{1,1,1,1},
- {1,0,0,1},
- {1,0,0,1},
- {1,1,1,1}};
- const unsigned char solid [4][4] = {{1,1,1,1},
- {1,1,1,1},
- {1,1,1,1},
- {1,1,1,1}};
- const unsigned char (*layers[4])[4][4] = {&solid, &hollow, &hollow, &solid};
-
- std::vector<CellData<dim> > cells;
- std::vector<bool> vertex_used (5*5*5, false);
-
- for (unsigned int zc = 0; zc < 4; ++zc)
- for (unsigned int yc = 0; yc < 4; ++yc)
- for (unsigned int xc = 0; xc < 4; ++xc)
+
+
+ // Radius of the sphere of the breast
+ // phantom geometry
+ static const double hemisphere_radius = 5;
+
+ // Radius of the disk underneath
+ static const double bottom_disk_radius = 10;
+
+ // Bottom z-coordinate of the disk
+ // underneath
+ const double bottom_disk_floor = -3;
+ // Top z-coordinate of the disk
+ // underneath
+ const double bottom_disk_ceil = -.5;
+
+ // radius of the inner set of cells of
+ // the sphere geometry
+ const double interior_hemisphere_radius
+ = hemisphere_radius/(1.+std::sqrt(2.0));
+
+ template <int dim>
+ class SphereBoundary : public HyperBallBoundary<dim>
+ {
+ public:
+ SphereBoundary ()
+ :
+ HyperBallBoundary<dim> (Point<dim>(), hemisphere_radius)
+ {}
+ };
+
+
+ template <int dim>
+ class CylinderBoundary : public StraightBoundary<dim>
+ {
+ public:
+ typedef
+ typename Triangulation<dim>::line_iterator
+ line_iterator;
+
+ typedef
+ typename Triangulation<dim>::quad_iterator
+ quad_iterator;
+
+ typedef
+ typename Triangulation<dim>::face_iterator
+ face_iterator;
+
+ /**
+ * Constructor.
+ */
+ CylinderBoundary (const double radius)
+ :
+ radius (radius)
+ {}
+
+
+ virtual Point<dim>
+ get_new_point_on_line (const line_iterator &line) const;
+
+ virtual Point<dim>
+ get_new_point_on_quad (const quad_iterator &quad) const;
+
+ virtual void
+ get_intermediate_points_on_line (const line_iterator &line,
+ std::vector<Point<dim> > &points) const;
+
+ virtual void
+ get_intermediate_points_on_quad (const quad_iterator &quad,
+ std::vector<Point<dim> > &points) const;
+
+ virtual void
+ get_normals_at_vertices (const face_iterator &face,
+ typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
+
+ private:
+ const double radius;
+
+ void
+ get_intermediate_points_between_points (const Point<dim> &p0, const Point<dim> &p1,
+ std::vector<Point<dim> > &points) const;
+ };
+
+
+ template <>
+ Point<3>
+ CylinderBoundary<3>::
+ get_new_point_on_line (const line_iterator &line) const
+ {
+ const Point<3> middle = StraightBoundary<3>::get_new_point_on_line (line);
+ // project to boundary
+ Point<3> p(middle[0], middle[1], 0);
+ p *= radius/std::sqrt(p.square());
+
+ return Point<3> (p[0], p[1], middle[2]);
+ }
+
+
+ template<>
+ Point<3>
+ CylinderBoundary<3>::
+ get_new_point_on_quad (const quad_iterator &quad) const
+ {
+ Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+
+ // project to boundary
+ Point<3> p(middle[0], middle[1], 0);
+ p *= radius/std::sqrt(p.square());
+
+ return Point<3> (p[0], p[1], middle[2]);
+ }
+
+
+ template <int dim>
+ void
+ CylinderBoundary<dim>::
+ get_intermediate_points_on_line (const 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);
+ }
+
+
+ template <int dim>
+ void
+ CylinderBoundary<dim>::
+ get_intermediate_points_between_points (const Point<dim> &,
+ const Point<dim> &,
+ std::vector<Point<dim> > &) const
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+
+ template <>
+ void
+ CylinderBoundary<3>::
+ get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &,
+ std::vector<Point<3> > &) const
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+
+ template <int dim>
+ void
+ CylinderBoundary<dim>::
+ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &,
+ typename Boundary<dim>::FaceVertexNormals &) const
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+
+ void
+ create_coarse_grid (Triangulation<3> &coarse_grid)
+ {
+ const unsigned int dim = 3;
+
+ std::vector<Point<dim> > vertices;
+ std::vector<CellData<dim> > cells;
+ SubCellData sub_cell_data;
+
+ const unsigned char
+ bottom_cylinder_boundary_id = 10,
+ middle_cylinder_boundary_id = 11,
+ spherical_boundary_id = 12,
+ all_other_boundary_id = 13,
+ straight_nondirichlet_boundary = 14;
+
+
+ // first build up the cells of the
+ // cylinder
+ {
+ // the vertices in each plane of
+ // the cylinder are located on
+ // three concentric rings of radii
+ // interior_hemisphere_radius,
+ // hemisphere_radius, and
+ // bottom_disk_radius,
+ // respectively. first generate
+ // these three rings
+ const Point<3> ring_points[8] = { Point<3>(-1,0,0),
+ Point<3>(-1,-1,0) / std::sqrt(2.),
+ Point<3>(0,-1,0),
+ Point<3>(+1,-1,0) / std::sqrt(2.),
+ Point<3>(+1,0,0),
+ Point<3>(+1,+1,0) / std::sqrt(2.),
+ Point<3>(0,+1,0),
+ Point<3>(-1,+1,0) / std::sqrt(2.) };
+
+ // first the point in the middle
+ // and the rest of those on the
+ // upper surface
+ vertices.push_back (Point<3>(0,0,bottom_disk_ceil));
+ for (unsigned int ring=0; ring<3; ++ring)
+ for (unsigned int i=0; i<8; ++i)
+ vertices.push_back (ring_points[i] * (ring == 0 ?
+ interior_hemisphere_radius :
+ (ring == 1 ? hemisphere_radius :
+ bottom_disk_radius))
+ +
+ Point<3>(0,0,bottom_disk_ceil));
+
+ // then points on lower surface
+ vertices.push_back (Point<3>(0,0,bottom_disk_floor));
+ for (unsigned int ring=0; ring<3; ++ring)
+ for (unsigned int i=0; i<8; ++i)
+ vertices.push_back (ring_points[i] * (ring == 0 ?
+ interior_hemisphere_radius :
+ (ring == 1 ?
+ hemisphere_radius :
+ bottom_disk_radius))
+ +
+ Point<3>(0,0,bottom_disk_floor));
+
+ const unsigned int n_vertices_per_surface = 25;
+ Assert (vertices.size() == n_vertices_per_surface*2,
+ ExcInternalError());
+
+ // next create cells from these
+ // vertices. only store the
+ // vertices of the upper surface,
+ // the lower ones are the same
+ // +12
+ {
+ const unsigned int connectivity[20][4]
+ = { { 1, 2, 3, 0 }, // four cells in the center
+ { 3, 4, 5, 0 },
+ { 0, 5, 6, 7 },
+ { 1, 0, 7, 8 },
+
+ { 9, 10, 2, 1 }, // eight cells of inner ring
+ { 10, 11, 3, 2 },
+ { 11, 12, 4, 3 },
+ { 4, 12, 13, 5 },
+ { 5, 13, 14, 6 },
+ { 6, 14, 15, 7 },
+ { 8, 7, 15, 16 },
+ { 9, 1, 8, 16 },
+
+ { 17, 18, 10, 9 }, // eight cells of outer ring
+ { 18, 19, 11, 10 },
+ { 19, 20, 12, 11 },
+ { 12, 20, 21, 13 },
+ { 13, 21, 22, 14 },
+ { 14, 22, 23, 15 },
+ { 16, 15, 23, 24 },
+ { 17, 9, 16, 24 } };
+
+ // now create cells out of this
+ for (unsigned int i=0; i<20; ++i)
+ {
+ CellData<3> cell;
+ for (unsigned int j=0; j<4; ++j)
+ {
+ cell.vertices[j] = connectivity[i][j];
+ cell.vertices[j+4] = connectivity[i][j]+n_vertices_per_surface;
+ }
+ cell.material_id = 0;
+ cells.push_back (cell);
+ }
+ }
+
+ // associate edges and faces on the
+ // outer boundary with boundary
+ // indicator of the cylinder
+ // boundary indicator. do this the
+ // same way as above, just this
+ // time with faces (edges follow
+ // from this immediately. some
+ // edges are duplicated since they
+ // belong to more than one cell,
+ // but that doesn't harm us here)
+ {
+ const unsigned int connectivity[8][2]
+ = { { 17,18 }, { 18, 19 }, { 19, 20 }, { 20, 21 },
+ { 21,22 }, { 22, 23 }, { 23, 24 }, { 24, 17 }};
+
+ for (unsigned int i=0; i<8; ++i)
+ {
+ const CellData<2> face =
+ { { connectivity[i][0]+n_vertices_per_surface,
+ connectivity[i][1]+n_vertices_per_surface,
+ connectivity[i][1],
+ connectivity[i][0] },
+ bottom_cylinder_boundary_id };
+ sub_cell_data.boundary_quads.push_back (face);
+
+ const CellData<1> edges[4] =
+ { { { connectivity[i][0], connectivity[i][1] },
+ bottom_cylinder_boundary_id },
+ { { connectivity[i][0]+n_vertices_per_surface,
+ connectivity[i][1]+n_vertices_per_surface },
+ bottom_cylinder_boundary_id },
+ { { connectivity[i][0]+n_vertices_per_surface,
+ connectivity[i][0] },
+ bottom_cylinder_boundary_id },
+ { { connectivity[i][1]+n_vertices_per_surface,
+ connectivity[i][1] },
+ bottom_cylinder_boundary_id } };
+ for (unsigned int i=0; i<4; ++i)
+ sub_cell_data.boundary_lines.push_back (edges[i]);
+ }
+ }
+ }
+
+ // next build up the middle ring. for
+ // this, copy the first 17 vertices
+ // up to z=0
+ {
+ const unsigned int first_upper_vertex = vertices.size();
+
+ for (unsigned int i=0; i<17; ++i)
+ vertices.push_back (Point<3>(vertices[i][0], vertices[i][1], 0));
+
+ // next create cells from these
+ // vertices. only store the
+ // vertices of the lower surface,
+ // the lower ones are the same
+ // +first_upper_vertex
+ const unsigned int connectivity[12][4]
+ = { { 1, 2, 3, 0 }, // four cells in the center
+ { 3, 4, 5, 0 },
+ { 0, 5, 6, 7 },
+ { 1, 0, 7, 8 },
+
+ { 9, 10, 2, 1 }, // eight cells of ring
+ { 10, 11, 3, 2 },
+ { 11, 12, 4, 3 },
+ { 4, 12, 13, 5 },
+ { 5, 13, 14, 6 },
+ { 6, 14, 15, 7 },
+ { 8, 7, 15, 16 },
+ { 9, 1, 8, 16 }};
+ // now create cells out of this
+ for (unsigned int i=0; i<12; ++i)
{
- // Check if we have to create a cell
- if ((*layers[zc])[xc][yc] == 1)
+ CellData<3> cell;
+ for (unsigned int j=0; j<4; ++j)
{
- const unsigned int z_vert = 25;
- const unsigned int y_vert = 5;
- unsigned int zoffs = zc * z_vert;
- unsigned int yoffs = yc * y_vert;
- unsigned int base_vert = zoffs + yoffs + xc;
-
- CellData<dim> cell;
- cell.vertices[0] = base_vert;
- cell.vertices[1] = cell.vertices[0] + 1;
- cell.vertices[2] = cell.vertices[0] + y_vert;
- cell.vertices[3] = cell.vertices[1] + y_vert;
- cell.vertices[4] = cell.vertices[0] + z_vert;
- cell.vertices[5] = cell.vertices[1] + z_vert;
- cell.vertices[6] = cell.vertices[2] + z_vert;
- cell.vertices[7] = cell.vertices[3] + z_vert;
- cell.material_id = 0;
- cells.push_back (cell);
-
- // Now add entries to the list of used
- // vertices.
- for (unsigned int i = 0; i < 8; ++i)
- vertex_used[cell.vertices[i]] = true;
+ cell.vertices[j] = connectivity[i][j]+first_upper_vertex;
+ cell.vertices[j+4] = connectivity[i][j];
}
+ cell.material_id = 0;
+ cells.push_back (cell);
}
- // Now create vertices and renumber stuff;
- std::vector<Point<dim> > vertices;
- std::vector<unsigned int> vert_renumber (5*5*5, 0);
- const double scale = 0.5;
- unsigned int v_indx = 0;
-
- for (int zv = 0; zv < 5; ++zv)
- for (int yv = 0; yv < 5; ++yv)
- for (int xv = 0; xv < 5; ++xv)
+ // mark the 8 vertical edges with
+ // the correct boundary indicator
+ for (unsigned int i=0; i<8; ++i)
{
- Point<dim> p_new ((double)(xv-2) * scale,
- (double)(yv-2) * scale,
- (double)(zv-2) * scale);
-
- if (vertex_used[v_indx])
- {
- vert_renumber[v_indx] = vertices.size ();
- vertices.push_back (p_new);
- }
- v_indx++;
+ const CellData<1> edge = { { 9, 9+first_upper_vertex },
+ middle_cylinder_boundary_id };
+ sub_cell_data.boundary_lines.push_back (edge);
+ }
+ // likewise with the 8 tangential
+ // edges on the lower disk. the
+ // edges at the interface between
+ // the middle disk and the
+ // hemisphere are handled by the
+ // hemisphere boundary
+ for (unsigned int i=0; i<8; ++i)
+ {
+ const CellData<1> edge = { { 9+i, 9+(i+1)%8},
+ middle_cylinder_boundary_id };
+ sub_cell_data.boundary_lines.push_back (edge);
}
- // Finally renumber the vertex indices in the cells
- std::vector<CellData<dim> >::iterator cell_iterator;
- for (cell_iterator = cells.begin (); cell_iterator != cells.end ();
- ++cell_iterator)
+ // then assign face indicators
+ for (unsigned int i=0; i<8; ++i)
+ {
+ const CellData<2> face = { { 9+i,
+ 9+(i+1)%8,
+ 9+(i+1)%8+first_upper_vertex,
+ 9+i+first_upper_vertex},
+ middle_cylinder_boundary_id };
+ sub_cell_data.boundary_quads.push_back (face);
+ }
+ }
+
+ // the final part is setting the
+ // half-sphere on top of this
{
- for (unsigned int i = 0; i < 8; ++i)
- cell_iterator->vertices[i] =
- vert_renumber[cell_iterator->vertices[i]];
+ // add four cubes to the top of
+ // the inner four cells, as well
+ // as 8 to their outside
+ {
+ // mirror the first nine vertices
+ // above the surface, and scale
+ // them to a certain distance
+ // outward
+ const double rx = hemisphere_radius / (1+std::sqrt(3.0));
+ for (unsigned int i=0; i<9; ++i)
+ {
+ const Point<3> p (vertices[i][0],
+ vertices[i][1],
+ i == 0 ?
+ 1
+ :
+ std::max(std::fabs(vertices[i][0]),
+ std::fabs(vertices[i][1])));
+ vertices.push_back (p / std::sqrt(p.square()) * rx);
+ }
+ Assert (vertices.size() == 76, ExcInternalError());
+
+ // same with the next ring of
+ // vertices, except that they
+ // go to hemisphere_radius
+ for (unsigned int i=9; i<17; ++i)
+ {
+ Point<3> p (vertices[i][0],
+ vertices[i][1],
+ std::max(std::fabs(vertices[i][0]),
+ std::fabs(vertices[i][1])));
+ vertices.push_back (p / std::sqrt(p.square()) *
+ hemisphere_radius);
+ }
+ Assert (vertices.size() == 84, ExcInternalError());
+
+ // make 12 cells out of this
+ const unsigned int connectivity[12][4]
+ = { { 1, 2, 3, 0 }, // four cells in the center
+ { 3, 4, 5, 0 },
+ { 0, 5, 6, 7 },
+ { 1, 0, 7, 8 },
+
+ { 9, 10, 2, 1 }, // eight cells of inner ring
+ { 10, 11, 3, 2 },
+ { 11, 12, 4, 3 },
+ { 4, 12, 13, 5 },
+ { 5, 13, 14, 6 },
+ { 6, 14, 15, 7 },
+ { 8, 7, 15, 16 },
+ { 9, 1, 8, 16 },
+ };
+
+ for (unsigned int i=0; i<12; ++i)
+ {
+ CellData<3> cell;
+ for (unsigned int j=0; j<4; ++j)
+ {
+ cell.vertices[j] = connectivity[i][j]+67;
+ cell.vertices[j+4] = connectivity[i][j]+50;
+ }
+ cell.material_id = 0;
+ cells.push_back (cell);
+ }
+ }
+
+ // assign boundary indicators to
+ // the faces and edges of these
+ // cells
+ {
+ // these are the numbers of the
+ // vertices on the top surface
+ // of the cylinder, with one
+ // "wrap-around":
+ const unsigned int vertices[9] =
+ { 9, 10, 11, 12, 13, 14, 15, 16, 9 };
+ // their counter-parts are the
+ // same +67
+ for (unsigned int i=0; i<8; ++i)
+ {
+ // generate a face
+ const CellData<2> face =
+ { { vertices[i]+50, vertices[i+1]+50,
+ vertices[i+1]+67, vertices[i]+67 },
+ spherical_boundary_id };
+ sub_cell_data.boundary_quads.push_back (face);
+
+ // same for the faces
+ const CellData<1> edges[4] =
+ { { { vertices[i]+50, vertices[i+1]+50 },
+ spherical_boundary_id },
+ { { vertices[i]+67, vertices[i+1]+67 },
+ spherical_boundary_id },
+ { { vertices[i]+50, vertices[i]+67 },
+ spherical_boundary_id },
+ { { vertices[i+1]+50, vertices[i+1]+67 },
+ spherical_boundary_id } };
+ for (unsigned int j=0; j<4; ++j)
+ sub_cell_data.boundary_lines.push_back (edges[j]);
+ }
+ }
+
+
+ // finally top the building
+ // with four closing cells and
+ // the vertex at the top
+ {
+ vertices.push_back (Point<3> (0,0,hemisphere_radius));
+
+ const unsigned int connectivity[4][8]
+ = { { 59, 60, 61, 67, 51, 52, 53, 50 },
+ { 61, 62, 63, 67, 53, 54, 55, 50 },
+ { 67, 63, 64, 65, 50, 55, 56, 57 },
+ { 59, 67, 65, 66, 51, 50, 57, 58 }};
+
+ for (unsigned int i=0; i<4; ++i)
+ {
+ CellData<3> cell;
+ for (unsigned int j=0; j<8; ++j)
+ cell.vertices[j] = connectivity[i][j]+17;
+ cell.material_id = 0;
+ cells.push_back (cell);
+ }
+
+ // generate boundary
+ // information for these cells,
+ // too
+ for (unsigned int i=0; i<4; ++i)
+ {
+ const CellData<2> face =
+ { { connectivity[i][0]+17, connectivity[i][1]+17,
+ connectivity[i][2]+17, connectivity[i][3]+17 },
+ spherical_boundary_id };
+ sub_cell_data.boundary_quads.push_back (face);
+
+ const CellData<1> edges[4] =
+ { { { connectivity[i][0]+17, connectivity[i][1]+17 },
+ spherical_boundary_id },
+ { { connectivity[i][1]+17, connectivity[i][2]+17 },
+ spherical_boundary_id },
+ { { connectivity[i][2]+17, connectivity[i][3]+17 },
+ spherical_boundary_id },
+ { { connectivity[i][3]+17, connectivity[i][0]+17 },
+ spherical_boundary_id } };
+ for (unsigned int j=0; j<4; ++j)
+ sub_cell_data.boundary_lines.push_back (edges[j]);
+ }
+ }
}
+
- // Now create triangulation
- triangulation.create_triangulation (vertices,
- cells,
- SubCellData());
+ // finally generate a triangulation
+ // out of this
+ GridReordering<3>::reorder_cells (cells);
+ coarse_grid.create_triangulation_compatibility (vertices, cells,
+ sub_cell_data);
+
+ // then associate boundary objects
+ // with the different boundary
+ // indicators
+ static const CylinderBoundary<3>
+ bottom_cylinder_boundary (bottom_disk_radius);
+ static const CylinderBoundary<3>
+ middle_cylinder_boundary (hemisphere_radius);
+ static const SphereBoundary<3> sphere_boundary;
+
+ coarse_grid.set_boundary (bottom_cylinder_boundary_id,
+ bottom_cylinder_boundary);
+ coarse_grid.set_boundary (middle_cylinder_boundary_id,
+ middle_cylinder_boundary);
+ coarse_grid.set_boundary (spherical_boundary_id,
+ sphere_boundary);
+
+ for (Triangulation<dim>::active_cell_iterator cell=coarse_grid.begin_active();
+ cell != coarse_grid.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if ((cell->face(f)->boundary_indicator() == 0)
+ &&
+ (cell->face(f)->center()[2] >= (bottom_disk_floor+bottom_disk_ceil)/2))
+ cell->face(f)->set_boundary_indicator(straight_nondirichlet_boundary);
+ }
+}
+
+
+
+template <>
+void LaplaceProblem<3>::create_coarse_grid ()
+{
+ BreastPhantom::create_coarse_grid (triangulation);
}
{
deallog.depth_console (0);
- LaplaceProblem<2> laplace_problem;
+ LaplaceProblem<3> laplace_problem;
laplace_problem.run ();
}
catch (std::exception &exc)