From: kayser-herold Date: Sun, 8 Apr 2007 19:05:50 +0000 (+0000) Subject: Introduced a constant for the maximum polynomial degree X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b378b8270be576d7de9e9e2f4d75b1941a74f786;p=dealii-svn.git Introduced a constant for the maximum polynomial degree that is used in the computation. Furthermore a method for creating a hollow cube in 3D was written. A hollow cube should be a more interesting case for hp-adaptivity as it introduces some additional corner singularities. Furthermore it is the 3D equivalent to the hollow square used in 2D. git-svn-id: https://svn.dealii.org/trunk@14621 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 86c632da27..8c4bd15c2f 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -90,6 +90,8 @@ class LaplaceProblem Vector system_rhs; Timer distr, condense, hang, assemble, solver; + + const unsigned int max_degree; }; @@ -121,9 +123,10 @@ RightHandSide::value (const Point &p, template LaplaceProblem::LaplaceProblem () : - dof_handler (triangulation) + dof_handler (triangulation), + max_degree (dim == 2 ? 7 : 4) { - for (unsigned int degree=2; degree<(dim == 2 ? 8 : 5); ++degree) + for (unsigned int degree=2; degree<=max_degree; ++degree) { fe_collection.push_back (FE_Q(degree)); quadrature_collection.push_back (QGauss(degree+2)); @@ -180,10 +183,9 @@ void LaplaceProblem::setup_system () condense.reset(); condense.start(); hanging_node_constraints.condense (csp); - condense.stop(); - + condense.stop(); sparsity_pattern.copy_from (csp); - } + } system_matrix.reinit (sparsity_pattern); } @@ -293,7 +295,7 @@ void LaplaceProblem:: estimate_smoothness (Vector &smoothness_indicators) const { - const unsigned int N = (dim == 2 ? 7 : 4); + const unsigned int N = max_degree; // form all the Fourier vectors // that we want to @@ -644,8 +646,96 @@ void LaplaceProblem<2>::create_coarse_grid () template <> void LaplaceProblem<3>::create_coarse_grid () { - GridGenerator::hyper_cube (triangulation); - triangulation.refine_global (1); + 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 > cells; + std::vector 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) + { + // Check if we have to create a cell + if ((*layers[zc])[xc][yc] == 1) + { + 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 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; + } + } + + // Now create vertices and renumber stuff; + std::vector > vertices; + std::vector 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) + { + Point 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++; + } + + // Finally renumber the vertex indices in the cells + std::vector >::iterator cell_iterator; + for (cell_iterator = cells.begin (); cell_iterator != cells.end (); + ++cell_iterator) + { + for (unsigned int i = 0; i < 8; ++i) + cell_iterator->vertices[i] = + vert_renumber[cell_iterator->vertices[i]]; + } + + // Now create triangulation + triangulation.create_triangulation (vertices, + cells, + SubCellData()); }