Hi Peter,
Great to hear from you, and thank you so much for pointing this out! I
revisited the implementation and corrected the setup for ghost cells as you
suggested. Now, the DoFs between subdomains are handled correctly, and the
total number of DoFs remains consistent. Your observation was spot
on—thanks for catching that!
One question: Can I set boundary_id as usual or do we have add in
construction_data?
I am attaching here boundary_ids implementation along with
generate_mesh_with_matrix.
template <int dim>
void LaplaceProblem<dim>::set_boundary_ids(int Nx, int Ny, int Nz)
{
// Boundary IDs
const types::boundary_id inlet_boundary_id = 0;
const types::boundary_id outlet_boundary_id = 1;
// Iterate over all active cells
for (const auto &cell : triangulation_pft.active_cell_iterators())
{
// Only modify locally owned cells
if (cell->is_locally_owned())
{
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++
face)
{
// Check if the face is at the inlet (x = 0)
if (std::abs(cell->face(face)->center()[0]) < 1e-12) // Near x = 0
{
cell->face(face)->set_boundary_id(inlet_boundary_id);
}
// Check if the face is at the outlet (x = Nx)
else if (std::abs(cell->face(face)->center()[0] - Nx) < 1e-12) // Near x =
Nx
{
cell->face(face)->set_boundary_id(outlet_boundary_id);
}
}
}
}
}
/////This one with GHost cells is important
template <int dim>
void LaplaceProblem<dim>::generate_mesh_from_matrix(int Nx, int Ny, int Nz)
{
// Prepare the construction data object
TriangulationDescription::Description<dim, dim> construction_data;
construction_data.comm = _mpi_communicator;
construction_data.cell_infos.resize(1); // Level 0 cell info
// Define the number of subdivisions and domain size
std::array<unsigned int, dim> n_subdivision = {{static_cast<unsigned int>
(Nx),
static_cast<unsigned int>(Ny),
static_cast<unsigned int>(Nz)}};
std::array<double, dim> sizes = {{1.0, 1.0, 1.0}}; // Assume a unit cube
domain for simplicity
// Partition the mesh along the x-axis
const unsigned int mpi_rank = Utilities::MPI::this_mpi_process
(_mpi_communicator);
const unsigned int mpi_size = Utilities::MPI::n_mpi_processes
(_mpi_communicator);
const auto partitioning_x = Utilities::
create_evenly_distributed_partitioning(mpi_rank, mpi_size, n_subdivision[0
]);
if (!partitioning_x.is_empty())
{
const int start_x_mpi = partitioning_x.nth_index_in_set(0);
const int end_x_mpi = partitioning_x.nth_index_in_set(partitioning_x.
n_elements() - 1);
// Generate vertices
std::vector<Point<dim>> vertices;
std::map<std::tuple<int, int, int>, unsigned int> local_vertex_map;
unsigned int vertex_counter = 0;
for (int i = std::max(0, start_x_mpi - 1); i <= std::min<int>(n_subdivision[
0], end_x_mpi + 1); ++i)
for (int j = 0; j <= n_subdivision[1]; ++j)
for (int k = 0; k <= n_subdivision[2]; ++k)
{
Point<dim> p(sizes[0] / n_subdivision[0] * i,
sizes[1] / n_subdivision[1] * j,
sizes[2] / n_subdivision[2] * k);
vertices.push_back(p);
local_vertex_map[std::make_tuple(i, j, k)] = vertex_counter++;
}
construction_data.coarse_cell_vertices = vertices;
// Generate local cells and a single layer of ghost cells
for (int i = std::max(0, start_x_mpi - 1); i <= std::min<int>(n_subdivision[
0] - 1, end_x_mpi + 1); ++i)
for (int j = 0; j < n_subdivision[1]; ++j)
for (int k = 0; k < n_subdivision[2]; ++k)
{
int material_id = mat_matrix[j + k * n_subdivision[1]][i];
if (material_id == _mat_tau_1)
{
CellData<dim> coarse_cell;
coarse_cell.vertices[0] = local_vertex_map[std::make_tuple(i, j, k)];
coarse_cell.vertices[1] = local_vertex_map[std::make_tuple(i + 1, j, k)];
coarse_cell.vertices[2] = local_vertex_map[std::make_tuple(i, j + 1, k)];
coarse_cell.vertices[3] = local_vertex_map[std::make_tuple(i + 1, j + 1,
k)];
if (dim == 3)
{
coarse_cell.vertices[4] = local_vertex_map[std::make_tuple(i, j, k + 1)];
coarse_cell.vertices[5] = local_vertex_map[std::make_tuple(i + 1, j, k + 1
)];
coarse_cell.vertices[6] = local_vertex_map[std::make_tuple(i, j + 1, k + 1
)];
coarse_cell.vertices[7] = local_vertex_map[std::make_tuple(i + 1, j + 1, k +
1)];
}
coarse_cell.material_id = material_id;
construction_data.coarse_cells.push_back(coarse_cell);
// Map the coarse cell index
unsigned int cell_id = i + j * n_subdivision[0] + k * n_subdivision[0] *
n_subdivision[1];
construction_data.coarse_cell_index_to_coarse_cell_id.push_back(cell_id);
// Fully distributed cell data
TriangulationDescription::CellData<dim> cell_data;
cell_data.id = CellId(cell_id, {}).template to_binary<dim>();
cell_data.subdomain_id = (i >= start_x_mpi && i <= end_x_mpi) ? mpi_rank :
(i < start_x_mpi ? mpi_rank - 1 : mpi_rank + 1);
cell_data.level_subdomain_id = 0;
construction_data.cell_infos[0].push_back(cell_data);
}
}
}
// Create the fully distributed triangulation
triangulation_pft.create_triangulation(construction_data);
}
However, I’ve run into a new issue with the solution. Somehow, I’m getting
a zero solution despite following the approach in step 40. I suspect that
the system might not be assembling properly.
Do you have any insights or advice on working with parallel fully
distributed meshes? I would really appreciate your thoughts.
Best Regards
Rishabh
On Tue, Dec 10, 2024 at 9:35 PM Peter Munch <peterrmue...@gmail.com> wrote:
> Hi Rishabh,
>
> To me it looks like that you don't set up the ghost cells properly. See
> that I loop over local cells and a single layer of ghost cells:
> https://github.com/dealii/dealii/blob/b2108aabb2ba014b6d35029db19034bded6f135a/tests/fullydistributed_grids/create_manually_01.cc#L70-L74.
> Once this is fixed, dofs between subdomains will not be duplicated so that
> the number of dofs shouls stay constant.
>
> Say hi to Thomas and the others!
>
> Best,
> Peter
>
> On Tuesday, 10 December 2024 at 16:25:36 UTC+1 risha...@gmail.com wrote:
>
>> Dear Prof. Wolfgang,
>>
>> Thank you for your guidance regarding the use of
>> TriangulationDescription::CellData and for suggesting Peter Munch's work
>> on creating fully distributed triangulations. It was very insightful, and I
>> reviewed the example provided in the repository:
>>
>> https://github.com/dealii/dealii/blob/master/tests/fullydistributed_grids/create_manually_01.cc
>>
>>
>> Based on this reference, I implemented a method to generate a fully
>> distributed triangulation for a small microstructure with dimensions
>> 20x20x20.
>> template <int dim>
>> void LaplaceProblem<dim>::generate_mesh_from_matrix(int Nx, int Ny, int
>> Nz)
>> {
>> // Prepare the construction data object
>> TriangulationDescription::Description<dim, dim> construction_data;
>> construction_data.comm = _mpi_communicator;
>> construction_data.cell_infos.resize(1); // Level 0 cell info
>>
>> // Define the number of subdivisions and domain size
>> std::array<unsigned int, dim> n_subdivision = {{static_cast<unsigned int>
>> (Nx),
>> static_cast<unsigned int>(Ny),
>> static_cast<unsigned int>(Nz)}};
>>
>> //std::array<unsigned int, dim> n_subdivision = {{2,2,2}};
>> std::array<double, dim> sizes = {{1.0, 1.0, 1.0}}; // Assume a unit cube
>> domain for simplicity
>>
>> // Partition the mesh along the x-axis
>> const unsigned int mpi_rank = Utilities::MPI::this_mpi_process
>> (_mpi_communicator);
>> const unsigned int mpi_size = Utilities::MPI::n_mpi_processes
>> (_mpi_communicator);
>> const auto partitioning_x = Utilities::
>> create_evenly_distributed_partitioning(mpi_rank, mpi_size, n_subdivision[
>> 0]);
>>
>> if (!partitioning_x.is_empty())
>> {
>> const int start_x_mpi = partitioning_x.nth_index_in_set(0);
>> const int end_x_mpi = partitioning_x.nth_index_in_set(partitioning_x.
>> n_elements() - 1);
>> //print start_x and end_x
>> std::cout << "start_x: " << start_x << std::endl;
>> std::cout << "end_x: " << end_x << std::endl;
>>
>> // Generate vertices
>> std::vector<Point<dim>> vertices;
>> std::map<std::tuple<int, int, int>, unsigned int> local_vertex_map;
>> unsigned int vertex_counter = 0;
>>
>> for (int i = std::max(0, start_x_mpi - 1); i <= std::min<int>(
>> n_subdivision[0], end_x_mpi + 1); ++i)
>> for (int j = 0; j <= n_subdivision[1]; ++j)
>> for (int k = 0; k <= n_subdivision[2]; ++k)
>> {
>> Point<dim> p(sizes[0] / n_subdivision[0] * i,
>> sizes[1] / n_subdivision[1] * j,
>> sizes[2] / n_subdivision[2] * k);
>> vertices.push_back(p);
>> local_vertex_map[std::make_tuple(i, j, k)] = vertex_counter++;
>> }
>>
>> construction_data.coarse_cell_vertices = vertices;
>>
>> // Generate local cells
>> for (int i = start_x_mpi; i <= end_x_mpi; ++i)
>> for (int j = 0; j < n_subdivision[1]; ++j)
>> for (int k = 0; k < n_subdivision[2]; ++k)
>> {
>> int material_id = mat_matrix[j + k * n_subdivision[1]][i];
>> if (material_id == _mat_tau_1)
>> {
>> CellData<dim> coarse_cell;
>> coarse_cell.vertices[0] = local_vertex_map[std::make_tuple(i, j, k)];
>> coarse_cell.vertices[1] = local_vertex_map[std::make_tuple(i + 1, j, k)];
>> coarse_cell.vertices[2] = local_vertex_map[std::make_tuple(i, j + 1, k)];
>> coarse_cell.vertices[3] = local_vertex_map[std::make_tuple(i + 1, j + 1,
>> k)];
>> if (dim == 3)
>> {
>> coarse_cell.vertices[4] = local_vertex_map[std::make_tuple(i, j, k + 1)];
>> coarse_cell.vertices[5] = local_vertex_map[std::make_tuple(i + 1, j, k +
>> 1)];
>> coarse_cell.vertices[6] = local_vertex_map[std::make_tuple(i, j + 1, k +
>> 1)];
>> coarse_cell.vertices[7] = local_vertex_map[std::make_tuple(i + 1, j + 1,
>> k + 1)];
>> }
>>
>> coarse_cell.material_id = material_id;
>> construction_data.coarse_cells.push_back(coarse_cell);
>>
>> // Map the coarse cell index
>> unsigned int cell_id = i + j * n_subdivision[0] + k * n_subdivision[0] *
>> n_subdivision[1];
>> construction_data.coarse_cell_index_to_coarse_cell_id.push_back(cell_id);
>>
>> // Fully distributed cell data
>> TriangulationDescription::CellData<dim> cell_data;
>> cell_data.id = CellId(cell_id, {}).template to_binary<dim>();
>> cell_data.subdomain_id = mpi_rank;
>> cell_data.level_subdomain_id = 0;
>>
>> construction_data.cell_infos[0].push_back(cell_data);
>> }
>> }
>> }
>>
>> // Create the fully distributed triangulation
>> triangulation_pft.create_triangulation(construction_data);
>>
>> // DataOut to handle distributed meshes
>> DataOut<dim> data_out;
>> data_out.attach_triangulation(triangulation_pft);
>>
>> // Add dummy data to ensure mesh output is generated
>> Vector<double> dummy_solution(triangulation_pft.n_active_cells());
>> for (unsigned int i = 0; i < dummy_solution.size(); ++i)
>> dummy_solution[i] = 0.0;
>> data_out.add_data_vector(dummy_solution, "dummy");
>>
>> // Build the output
>> data_out.build_patches();
>>
>> // Output VTU file
>> std::string filename = "mesh-" + Utilities::int_to_string(mpi_rank, 4) +
>> ".vtu";
>> std::ofstream output(filename);
>> data_out.write_vtu(output);
>>
>> }
>>
>>
>>
>> The VTU files generated for each processor appear correct and are divided
>> as expected. However, I am encountering an issue during the DoF
>> distribution phase.
>> Specifically, the number of DoFs increases as I use more processors. For
>> instance, using a single processor results in 7918 DoFs, while with two
>> processors, the count rises to 8272, and with three processors, it
>> increases further to 8537.
>>
>> This suggests that there might be an overlapping of cell faces or an
>> issue with ghost cells being incorrectly accounted for during the DoF
>> distribution. I suspect the problem lies in the handling of shared faces
>> between processors or inconsistencies in cell ownership. I am currently
>> reviewing the partitioning logic and cell indexing within
>> TriangulationDescription::Description.
>>
>> If you have any further suggestions or recommendations to address this
>> issue, they would be greatly appreciated. Thank you once again for your
>> support, and I look forward to your feedback.
>>
>> Best regards,
>> Rishabh
>>
>> On Tue, Dec 10, 2024 at 2:03 PM Wolfgang Bangerth <bang...@colostate.edu>
>> wrote:
>>
>>> On 12/9/24 08:46, Rishabh Saxena wrote:
>>> > cell.vertices = {vertex_map[std::make_tuple(i, j,
>>> k)],
>>>
>>> >[...]
>>>
>>> > When I compile the code, I get the following errors:
>>> >
>>> > /home/saxenar/local/programs/tortuosity_github/source/
>>> > tortuosity_parallel_FDM_discriptor_MPI.cc:957:26: Fehler: »struct
>>> > dealii::TriangulationDescription::CellData<3>« hat kein Element namens
>>> »vertices«
>>>
>>> Rishabh:
>>> the error message says that TriangulationDescription::CellData simply
>>> does not
>>> have member variables 'vertices', 'material_id', etc., and that is true:
>>>
>>>
>>> https://dealii.org/developer/doxygen/deal.II/structTriangulationDescription_1_1CellData.html
>>> The general description of the class also explains why it does not need
>>> to
>>> store the geometric information you are trying to describe.
>>>
>>> I believe that Peter Munch posted a piece of code here on this forum a
>>> while
>>> ago that creates an xyz-uniform fully distributed triangulation like the
>>> one
>>> you are setting up here. You might want to search through the archive
>>> and see
>>> if you can find it!
>>>
>>> Best
>>> WB
>>>
>>> --
>>> 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 dealii+un...@googlegroups.com.
>>> To view this discussion visit
>>> https://groups.google.com/d/msgid/dealii/14fb678f-ffb8-40c1-abce-29bec8078fee%40colostate.edu
>>> .
>>>
>>
>>
>> --
>>
>>
>>
>>
>> *Thanks & RegardsRishabh SaxenaHelmut Schmidt University || University of
>> the Federal Armed Forces HamburgApplied Mathematics Holstenhofweg 85 22
>> 043 Hamburg GermanyOffice: H 1, Room No. 1225Email: sax...@hsu-hh.de*
>>
> --
> 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 dealii+unsubscr...@googlegroups.com.
> To view this discussion visit
> https://groups.google.com/d/msgid/dealii/c703ffab-89b5-4e7d-b363-10dc062ec9fan%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/c703ffab-89b5-4e7d-b363-10dc062ec9fan%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
--
*Thanks & RegardsRishabh SaxenaHelmut Schmidt University || University of
the Federal Armed Forces HamburgApplied Mathematics Holstenhofweg 85 22
043 Hamburg GermanyOffice: H 1, Room No. 1225Email: saxe...@hsu-hh.de
<saxe...@hsu-hh.de>*
--
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 dealii+unsubscr...@googlegroups.com.
To view this discussion visit
https://groups.google.com/d/msgid/dealii/CADCJbzORxFfbU3v1vgPxdAgQJHffi_OB9MrP0mVAw1%2BaTvg7DA%40mail.gmail.com.
