#include <deal.II/lac/trilinos_sparsity_pattern.h>
#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
#include <deal.II/simplex/fe_lib.h>
#include <deal.II/simplex/grid_generator.h>
SolverCG<VectorType> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+ // deallog << solution.l2_norm() << " " << system_rhs.l2_norm() << std::endl;
+
deallog << " with " << solver_control.last_step()
<< " CG iterations needed to obtain convergence" << std::endl;
deallog << std::endl;
+
if (false)
{
+ Vector<double> difference(tria.n_active_cells());
+
+ VectorTools::integrate_difference(mapping,
+ dof_handler,
+ solution,
+ Functions::ZeroFunction<dim>(),
+ difference,
+ quad,
+ VectorTools::NormType::L2_norm);
+
+ deallog << VectorTools::compute_global_error(
+ tria, difference, VectorTools::NormType::L2_norm)
+ << std::endl;
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
test(*tria, fe, quad, face_quad, mapping, params.p2[0], false);
}
+template <int dim, int spacedim = dim>
+void
+test_pyramid(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_pyramids(
+ *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_PyramidP<dim> fe(params.degree);
+
+ Simplex::QGaussPyramid<dim> quad(params.degree + 1);
+
+ Quadrature<dim - 1> face_quad; // not needed
+
+ Simplex::FE_PyramidP<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>(3.2, 1, 1);
test_wedge(comm, params);
}
+
+ // test PYRAMID
+ {
+ deallog << "Solve problem on PYRAMID mesh:" << std::endl;
+
+ params.file_name_out = "mesh-pyramid";
+ params.repetitions = std::vector<unsigned int>{10, 10, 10};
+ params.p1 = Point<3>(3.3, 0, 0);
+ params.p2 = Point<3>(4.3, 1, 1);
+
+ params.degree = 1;
+
+ test_pyramid(comm, params);
+ }
}
}
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_pyramids(
+ 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));
+ 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 + 0.5),
+ p1[1] + dx[1] * (j + 0.5),
+ p1[2] + dx[2] * (k + 0.5)));
+
+ // 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, 9> 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,
+ (repetitions[2] + 1) * (repetitions[1] + 1) *
+ (repetitions[0] + 1) +
+ (repetitions[1] * repetitions[0] * k +
+ repetitions[0] * j + i)}};
+
+
+ // TRI cell 0
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[0], quad[2], quad[4], quad[6], quad[8]};
+ cells.push_back(tri);
+ }
+ // TRI cell 1
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[2], quad[3], quad[6], quad[7], quad[8]};
+ cells.push_back(tri);
+ }
+ // TRI cell 2
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[6], quad[7], quad[4], quad[5], quad[8]};
+ cells.push_back(tri);
+ }
+ // TRI cell 3
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[1], quad[5], quad[3], quad[7], quad[8]};
+ cells.push_back(tri);
+ }
+ // TRI cell 4
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[0], quad[1], quad[2], quad[3], quad[8]};
+ cells.push_back(tri);
+ }
+ // TRI cell 5
+ {
+ CellData<dim> tri;
+ tri.vertices = {
+ quad[0], quad[4], quad[1], quad[5], quad[8]};
+ 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