#include "../tests.h"
+#include "simplex_grids.h"
+
using namespace dealii;
template <int dim>
const Quadrature<dim> & quad,
const Quadrature<dim - 1> & face_quad,
const Mapping<dim, spacedim> & mapping,
- const double r_boundary)
+ const double r_boundary,
+ const bool do_use_fe_face_values = true)
{
std::string label =
(dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
std::shared_ptr<FEFaceValues<dim, spacedim>> fe_face_values;
- fe_face_values.reset(
- new FEFaceValues<dim, spacedim>(mapping, fe, face_quad, flag));
+ if (do_use_fe_face_values)
+ fe_face_values.reset(
+ new FEFaceValues<dim, spacedim>(mapping, fe, face_quad, flag));
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quad.size();
hex_mesh &= (cell->n_vertices() == GeometryInfo<dim>::vertices_per_cell);
deallog << std::endl;
+
+ if (false)
+ {
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "solution");
+
+ data_out.build_patches(mapping);
+
+ std::ofstream output("result.vtk");
+ data_out.write_vtk(output);
+ }
}
template <int dim, int spacedim = dim>
test(tria, fe, quad, quad_face, mapping, params.p2[0]);
}
+template <int dim, int spacedim = dim>
+void
+test_wedge(const MPI_Comm &comm, const Parameters<dim> ¶ms)
+{
+ const unsigned int tria_type = 2;
+
+ // 1) Create triangulation...
+ Triangulation<dim, spacedim> *tria;
+
+ // a) serial triangulation
+ Triangulation<dim, spacedim> tr_1;
+
+ // b) shared triangulation (with artificial cells)
+ parallel::shared::Triangulation<dim> tr_2(
+ MPI_COMM_WORLD,
+ ::Triangulation<dim>::none,
+ true,
+ parallel::shared::Triangulation<dim>::partition_custom_signal);
+
+ tr_2.signals.create.connect([&]() {
+ GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm),
+ tr_2);
+ });
+
+ // c) distributed triangulation
+ parallel::fullydistributed::Triangulation<dim> tr_3(comm);
+
+
+ // ... choose the right triangulation
+ if (tria_type == 0 || tria_type == 2)
+ tria = &tr_1;
+ else if (tria_type == 1)
+ tria = &tr_2;
+
+ // ... create triangulation
+ if (params.use_grid_generator)
+ {
+ // ...via Simplex::GridGenerator
+ GridGenerator::subdivided_hyper_rectangle_with_wedges(
+ *tria, params.repetitions, params.p1, params.p2, false);
+ }
+ else
+ {
+ // ...via GridIn
+ GridIn<dim, spacedim> grid_in;
+ grid_in.attach_triangulation(*tria);
+ std::ifstream input_file(params.file_name_in);
+ grid_in.read_ucd(input_file);
+ // std::ifstream input_file("test_tet_geometry.unv");
+ // grid_in.read_unv(input_file);
+ }
+
+ // ... partition serial triangulation and create distributed triangulation
+ if (tria_type == 0 || tria_type == 2)
+ {
+ GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm),
+ tr_1);
+
+ auto construction_data = TriangulationDescription::Utilities::
+ create_description_from_triangulation(tr_1, comm);
+
+ tr_3.create_triangulation(construction_data);
+
+ tria = &tr_3;
+ }
+
+ // 2) Output generated triangulation via GridOut
+ GridOut grid_out;
+ std::ofstream out(params.file_name_out + "." +
+ std::to_string(Utilities::MPI::this_mpi_process(comm)) +
+ ".vtk");
+ grid_out.write_vtk(*tria, out);
+
+ // 3) Select components
+ Simplex::FE_WedgeP<dim> fe(params.degree);
+
+ Simplex::QGaussWedge<dim> quad(params.degree + 1);
+
+ Quadrature<dim - 1> face_quad; // not needed
+
+ Simplex::FE_WedgeP<dim> fe_mapping(1);
+ MappingFE<dim> mapping(fe_mapping);
+
+ // 4) Perform test (independent of mesh type)
+ test(*tria, fe, quad, face_quad, mapping, params.p2[0], false);
+}
+
int
main(int argc, char **argv)
{
params.p2 = Point<3>(2.1, 1, 1);
test_hex(comm, params);
}
+
+ // test WEDGE
+ {
+ deallog << "Solve problem on WEDGE mesh:" << std::endl;
+
+ params.file_name_out = "mesh-wedge";
+ params.p1 = Point<3>(2.2, 0, 0);
+ params.p2 = Point<3>(3.2, 1, 1);
+ test_wedge(comm, params);
+ }
}
}
DEAL::Solve problem on TRI mesh:
DEAL:: on parallel::fullydistributed::Triangulation
DEAL:cg::Starting value 0.245798
-DEAL:cg::Convergence step 114 value 6.95996e-13
+DEAL:cg::Convergence step 114 value 6.62438e-13
DEAL:: with 114 CG iterations needed to obtain convergence
DEAL::
DEAL::Solve problem on QUAD mesh:
DEAL::Solve problem on TET mesh:
DEAL:: on parallel::fullydistributed::Triangulation
DEAL:cg::Starting value 0.0607616
-DEAL:cg::Convergence step 156 value 9.07635e-13
+DEAL:cg::Convergence step 156 value 9.07675e-13
DEAL:: with 156 CG iterations needed to obtain convergence
DEAL::
DEAL::Solve problem on HEX mesh:
DEAL:cg::Convergence step 134 value 8.23917e-13
DEAL:: with 134 CG iterations needed to obtain convergence
DEAL::
+DEAL::Solve problem on WEDGE mesh:
+DEAL:: on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.0132550
+DEAL:cg::Convergence step 186 value 8.72381e-13
+DEAL:: with 186 CG iterations needed to obtain convergence
+DEAL::
DEAL::Solve problem on TRI mesh:
DEAL:: on parallel::fullydistributed::Triangulation
DEAL:cg::Starting value 0.245798
-DEAL:cg::Convergence step 114 value 6.45793e-13
+DEAL:cg::Convergence step 114 value 6.64515e-13
DEAL:: with 114 CG iterations needed to obtain convergence
DEAL::
DEAL::Solve problem on QUAD mesh:
DEAL::Solve problem on TET mesh:
DEAL:: on parallel::fullydistributed::Triangulation
DEAL:cg::Starting value 0.0607616
-DEAL:cg::Convergence step 156 value 9.07640e-13
+DEAL:cg::Convergence step 156 value 9.07649e-13
DEAL:: with 156 CG iterations needed to obtain convergence
DEAL::
DEAL::Solve problem on HEX mesh:
DEAL:cg::Convergence step 134 value 8.23917e-13
DEAL:: with 134 CG iterations needed to obtain convergence
DEAL::
+DEAL::Solve problem on WEDGE mesh:
+DEAL:: on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.0132550
+DEAL:cg::Convergence step 186 value 8.72414e-13
+DEAL:: with 186 CG iterations needed to obtain convergence
+DEAL::
{
namespace GridGenerator
{
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_wedges(
+ Triangulation<dim, spacedim> & tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<dim> & p1,
+ const Point<dim> & p2,
+ const bool colorize = false)
+ {
+ AssertDimension(dim, spacedim);
+
+ AssertThrow(colorize == false, ExcNotImplemented());
+
+ std::vector<Point<spacedim>> vertices;
+ std::vector<CellData<dim>> cells;
+
+ if (dim == 3)
+ {
+ // determine cell sizes
+ const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+ (p2[1] - p1[1]) / repetitions[1],
+ (p2[2] - p1[2]) / repetitions[2]);
+
+ // create vertices
+ for (unsigned int k = 0; k <= repetitions[2]; ++k)
+ for (unsigned int j = 0; j <= repetitions[1]; ++j)
+ for (unsigned int i = 0; i <= repetitions[0]; ++i)
+ vertices.push_back(Point<spacedim>(p1[0] + dx[0] * i,
+ p1[1] + dx[1] * j,
+ p1[2] + dx[2] * k));
+
+ // create cells
+ for (unsigned int k = 0; k < repetitions[2]; ++k)
+ for (unsigned int j = 0; j < repetitions[1]; ++j)
+ for (unsigned int i = 0; i < repetitions[0]; ++i)
+ {
+ // create reference HEX cell
+ std::array<unsigned int, 8> quad{
+ {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 0,
+ (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 1,
+ (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 0,
+ (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 1,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 0,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 1,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 0,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 1}};
+
+
+ // TRI cell 0
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[0], quad[1], quad[2], quad[4], quad[5], quad[6]};
+ cells.push_back(tri);
+ }
+
+ // TRI cell 1
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[3], quad[2], quad[1], quad[7], quad[6], quad[5]};
+ cells.push_back(tri);
+ }
+ }
+ }
+ else
+ {
+ AssertThrow(colorize == false, ExcNotImplemented());
+ }
+
+ // actually create triangulation
+ tria.create_triangulation(vertices, cells, SubCellData());
+ }
+
template <int dim, int spacedim>
void
subdivided_hyper_rectangle_with_simplices_mix(
projects / dealii.git / commitdiff