# Mixed topology

Experimental demo.

#include <dolfinx/common/IndexMap.h>
#include <dolfinx/common/MPI.h>
#include <dolfinx/mesh/Geometry.h>
#include <dolfinx/mesh/Mesh.h>
#include <dolfinx/mesh/Topology.h>
#include <dolfinx/mesh/cell_types.h>
#include <iostream>
#include <vector>

using namespace dolfinx;


Note: this demo is not currently intended to provide a fully functional example of using a mixed-topology mesh, but shows only the basic construction. Experimental.

int main(int argc, char* argv[])
{
dolfinx::init_logging(argc, argv);
MPI_Init(&argc, &argv);

// Number of square cell in x-direction
constexpr int nx_s = 2;
// Number of triangle cells in x-direction
constexpr int nx_t = 2;
// Number of cells in y-direction
constexpr int ny = 4;

constexpr int num_s = nx_s * ny;
constexpr int num_t = 2 * nx_t * ny;

std::vector<double> x;
for (int i = 0; i < nx_s + nx_t + 1; ++i)
{
for (int j = 0; j < ny + 1; ++j)
{
x.push_back(i);
x.push_back(j);
}
}

std::vector<std::int64_t> cells;
std::vector<std::int32_t> offsets{0};
for (int i = 0; i < nx_s + nx_t; ++i)
{
for (int j = 0; j < ny; ++j)
{
const int v_0 = j + i * (ny + 1);
const int v_1 = v_0 + 1;
const int v_2 = v_0 + ny + 1;
const int v_3 = v_0 + ny + 2;
if (i < nx_s)
{
cells.push_back(v_0);
cells.push_back(v_1);
cells.push_back(v_3);
cells.push_back(v_2);
offsets.push_back(cells.size());
}
else
{
cells.push_back(v_0);
cells.push_back(v_1);
cells.push_back(v_2);
offsets.push_back(cells.size());

cells.push_back(v_1);
cells.push_back(v_2);
cells.push_back(v_3);
offsets.push_back(cells.size());
}
}
}

std::vector<std::int64_t> original_global_index(num_s + num_t);
std::iota(original_global_index.begin(), original_global_index.end(), 0);
std::vector<int> ghost_owners;
std::vector<std::int32_t> cell_group_offsets{0, num_s, num_s + num_t,
num_s + num_t, num_s + num_t};
std::vector<std::int64_t> boundary_vertices;
for (int j = 0; j < ny + 1; ++j)
{
boundary_vertices.push_back(j);
boundary_vertices.push_back(j + (nx_s + nx_t + 1) * ny);
}
for (int i = 0; i < nx_s + nx_t + 1; ++i)
{
boundary_vertices.push_back((ny + 1) * i);
boundary_vertices.push_back(ny + (ny + 1) * i);
}

std::sort(boundary_vertices.begin(), boundary_vertices.end());
boundary_vertices.erase(
std::unique(boundary_vertices.begin(), boundary_vertices.end()),
boundary_vertices.end());

mesh::CellType::triangle};
std::vector<fem::CoordinateElement<double>> elements;
for (auto ct : cell_types)
elements.push_back(fem::CoordinateElement<double>(ct, 1));

{
auto topo = std::make_shared<mesh::Topology>(mesh::create_topology(
MPI_COMM_WORLD, cells_list, original_global_index, ghost_owners,
cell_types, cell_group_offsets, boundary_vertices));

auto topo_cells = topo->connectivity(2, 0);

std::cout << "cells\n------\n";
for (int i = 0; i < topo_cells->num_nodes(); ++i)
{
std::cout << i << " [";
std::cout << q << " ";
std::cout << "]\n";
}

topo->create_connectivity(1, 0);

std::cout << "facets\n------\n";
auto topo_facets = topo->connectivity(1, 0);
for (int i = 0; i < topo_facets->num_nodes(); ++i)
{
std::cout << i << " [";
std::cout << q << " ";
std::cout << "]\n";
}

mesh::Geometry geom = mesh::create_geometry(MPI_COMM_WORLD, *topo, elements,
cells_list, x, 2);

mesh::Mesh<double> mesh(MPI_COMM_WORLD, topo, geom);
std::cout << "num cells = " << mesh.topology()->index_map(2)->size_local()
<< "\n";
for (auto q : mesh.topology()->entity_group_offsets(2))
std::cout << q << " ";
std::cout << "\n";
}

MPI_Finalize();

return 0;
}