Dear all,
I faced a strange problem. With relatively few cells and non-conforming
triangulation, the assembling process generates new nonzero locations in
parallel with relatively large amount of processes . It works well in
serial or with a small number of processes or with conforming
triangulation. I wonder whether it is a bug or my fault. Can anybody
provide some hints?
The attached file is a minimal example to reproduce the problem. Please run
with 4 or more processes.
Best regards,
Ce Qin
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
int main(int argc, char **argv) {
dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
size_t n_mpi_processes = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
size_t this_mpi_process = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
dealii::FE_Q<3> fe(1);
dealii::Triangulation<3> tria;
dealii::GridGenerator::subdivided_hyper_cube(tria, 2);
tria.begin_active()->set_refine_flag();
tria.execute_coarsening_and_refinement();
dealii::GridTools::partition_triangulation(n_mpi_processes, tria);
dealii::DoFHandler<3> dh(tria);
dh.distribute_dofs(fe);
dealii::DoFRenumbering::subdomain_wise(dh);
if (this_mpi_process == 0) {
std::cout << "Number of Cells: " << tria.n_active_cells() << std::endl;
std::cout << "Number of DoFs: " << dh.n_dofs() << std::endl;
}
dealii::ConstraintMatrix cm;
dealii::DoFTools::make_hanging_node_constraints(dh, cm);
dealii::SparsityPattern sp;
dealii::DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs());
dealii::DoFTools::make_sparsity_pattern(dh, dsp);
sp.copy_from(dsp);
std::vector<dealii::types::global_dof_index> local_rows_per_process(n_mpi_processes);
for (size_t p = 0; p < n_mpi_processes; ++p) {
local_rows_per_process[p] = dealii::DoFTools::count_dofs_with_subdomain_association(dh, p);
if (this_mpi_process == 0) {
std::cout << p << " : " << local_rows_per_process[p] << std::endl;
}
}
dealii::PETScWrappers::MPI::SparseMatrix system_matrix;
system_matrix.reinit(MPI_COMM_WORLD, sp, local_rows_per_process, local_rows_per_process,
this_mpi_process);
dealii::QGauss<3> quadrature_formula(2);
dealii::FEValues<3> fe_values(fe, quadrature_formula, dealii::update_values |
dealii::update_quadrature_points |
dealii::update_JxW_values);
const size_t dofs_per_cell = fe.dofs_per_cell;
const size_t n_q_points = quadrature_formula.size();
dealii::FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
std::vector<dealii::types::global_dof_index> local_dof_indices(dofs_per_cell);
for (auto cell = dh.begin_active(); cell != dh.end(); ++cell) {
if (cell->subdomain_id() != this_mpi_process) {
continue;
}
fe_values.reinit(cell);
cell_matrix = 0.0;
for (size_t i = 0; i < dofs_per_cell; ++i) {
for (size_t j = 0; j < dofs_per_cell; ++j) {
for (size_t q_point = 0; q_point < n_q_points; ++q_point) {
cell_matrix(i, j) +=
(fe_values.shape_value(i, q_point) * fe_values.shape_value(j, q_point)) *
fe_values.JxW(q_point);
}
}
}
cell->get_dof_indices(local_dof_indices);
cm.distribute_local_to_global(cell_matrix, local_dof_indices, system_matrix);
}
system_matrix.compress(dealii::VectorOperation::add);
return 0;
}
