From 1c09ea65a9c6140cb8d017279d42483e06a44375 Mon Sep 17 00:00:00 2001 From: "lisa.pizzo" Date: Sat, 10 Jan 2026 15:26:59 +0100 Subject: [PATCH] . --- .../jacob_template/ascii_read_meshvector.m | 43 ++++ Sheet7/E14/jacob_template/ascii_write_mesh.m | 54 +++++ .../jacob_template/ascii_write_subdomains.m | 51 ++++ .../jacob_template/cuthill_mckee_ordering.cpp | 218 ++++++++++++++++++ .../jacob_template/cuthill_mckee_ordering.h | 8 + 5 files changed, 374 insertions(+) create mode 100644 Sheet7/E14/jacob_template/ascii_read_meshvector.m create mode 100644 Sheet7/E14/jacob_template/ascii_write_mesh.m create mode 100644 Sheet7/E14/jacob_template/ascii_write_subdomains.m create mode 100644 Sheet7/E14/jacob_template/cuthill_mckee_ordering.cpp create mode 100644 Sheet7/E14/jacob_template/cuthill_mckee_ordering.h diff --git a/Sheet7/E14/jacob_template/ascii_read_meshvector.m b/Sheet7/E14/jacob_template/ascii_read_meshvector.m new file mode 100644 index 0000000..03119db --- /dev/null +++ b/Sheet7/E14/jacob_template/ascii_read_meshvector.m @@ -0,0 +1,43 @@ +function [ xc, ia, v ] = ascii_read_meshvector( fname ) +% +% Loads the 2D triangular mesh (coordinates, vertex connectivity) +% together with values on its vertices from an ASCII file. +% Matlab indexing is stored (starts with 1). +% +% The input file format is compatible +% with Mesh_2d_3_matlab:Write_ascii_matlab(..) in jacobi_oo_stl/geom.h +% +% +% IN: fname - filename +% OUT: xc - coordinates +% ia - mesh connectivity +% v - solution vector + +DELIMETER = ' '; + +fprintf('Read file %s\n',fname) + +% Read mesh constants +nn = dlmread(fname,DELIMETER,[0 0 0 3]); %% row_1, col_1, row_2, col_2 in C indexing!!! +nnode = nn(1); +ndim = nn(2); +nelem = nn(3); +nvert = nn(4); + +% Read coordinates +row_start = 0+1; +row_end = 0+nnode; +xc = dlmread(fname,DELIMETER,[row_start 0 row_end ndim-1]); + +% Read connectivity +row_start = row_end+1; +row_end = row_end+nelem; +ia = dlmread(fname,DELIMETER,[row_start 0 row_end nvert-1]); + +% Read solution +row_start = row_end+1; +row_end = row_end+nnode; +v = dlmread(fname,DELIMETER,[row_start 0 row_end 0]); +end + + diff --git a/Sheet7/E14/jacob_template/ascii_write_mesh.m b/Sheet7/E14/jacob_template/ascii_write_mesh.m new file mode 100644 index 0000000..262fef8 --- /dev/null +++ b/Sheet7/E14/jacob_template/ascii_write_mesh.m @@ -0,0 +1,54 @@ +function ascii_write_mesh( xc, ia, e, basename) +% +% Saves the 2D triangular mesh in the minimal way (only coordinates, vertex connectivity, minimal boundary edge info) +% in an ASCII file. +% Matlab indexing is stored (starts with 1). +% +% The output file format is compatible with Mesh_2d_3_matlab:Mesh_2d_3_matlab(std::string const &fname) in jacobi_oo_stl/geom.h +% +% IN: +% coordinates xc: [2][nnode] +% connectivity ia: [4][nelem] with t(4,:) are the subdomain numbers +% edges e: [7][nedges] boundary edges +% e([1,2],:) - start/end vertex of edge +% e([3,4],:) - start/end values +% e(5,:) - segment number +% e([6,7],:) - left/right subdomain +% basename: file name without extension +% +% Data have been generated via . +% +fname = [basename, '.txt']; + +nnode = int32(size(xc,2)); +ndim = int32(size(xc,1)); +nelem = int32(size(ia,2)); +nvert_e = int32(3); + + +dlmwrite(fname,nnode,'delimiter','\t','precision',16) % number of nodes +dlmwrite(fname,ndim,'-append','delimiter','\t','precision',16) % space dimension +dlmwrite(fname,nelem,'-append','delimiter','\t','precision',16) % number of elements +dlmwrite(fname,nvert_e,'-append','delimiter','\t','precision',16) % number of vertices per element + +%% dlmwrite(fname,xc(:),'-append','delimiter','\t','precision',16) % coordinates +dlmwrite(fname,xc([1,2],:).','-append','delimiter','\t','precision',16) % coordinates + +%% no subdomain info transferred +tmp=int32(ia(1:3,:)); +% dlmwrite(fname,tmp(:),'-append','delimiter','\t','precision',16) % connectivity in Matlab indexing +dlmwrite(fname,tmp(:,:).','-append','delimiter','\t','precision',16) % connectivity in Matlab indexing + +%% store only start and end point of boundary edges, +nbedges = size(e,2); +dlmwrite(fname,nbedges,'-append','delimiter','\t','precision',16) % number boundary edges +tmp=int32(e(1:2,:)); +% dlmwrite(fname,tmp(:),'-append','delimiter','\t','precision',16) % boundary edges in Matlab indexing +dlmwrite(fname,tmp(:,:).','-append','delimiter','\t','precision',16) % boundary edges in Matlab indexing + +%% Add subdomain information to edges +tmp=int32(e(6:7,:)); +dlmwrite(fname,tmp(:,:).','-append','delimiter','\t','precision',16) % subdomain information to edges + + +end diff --git a/Sheet7/E14/jacob_template/ascii_write_subdomains.m b/Sheet7/E14/jacob_template/ascii_write_subdomains.m new file mode 100644 index 0000000..dd9c6a7 --- /dev/null +++ b/Sheet7/E14/jacob_template/ascii_write_subdomains.m @@ -0,0 +1,51 @@ +function ascii_write_subdomains( xc, ia, e, basename) +% +% Saves the 2D triangular mesh in the minimal way (only coordinates, vertex connectivity, minimal boundary edge info) +% in an ASCII file. +% Matlab indexing is stored (starts with 1). +% +% The output file format is compatible with Mesh_2d_3_matlab:Mesh_2d_3_matlab(std::string const &fname) in jacobi_oo_stl/geom.h +% +% IN: +% coordinates xc: [2][nnode] +% connectivity ia: [4][nelem] with t(4,:) are the subdomain numbers +% edges e: [7][nedges] boundary edges +% e([1,2],:) - start/end vertex of edge +% e([3,4],:) - start/end values +% e(5,:) - segment number +% e([6,7],:) - left/right subdomain +% basename: file name without extension +% +% Data have been generated via . +% +fname = [basename, '_sd.txt']; + +nnode = int32(size(xc,2)); +ndim = int32(size(xc,1)); +nelem = int32(size(ia,2)) +nvert_e = int32(3); + + +% dlmwrite(fname,nnode,'delimiter','\t','precision',16) % number of nodes +% dlmwrite(fname,ndim,'-append','delimiter','\t','precision',16) % space dimension +% dlmwrite(fname,nelem,'-append','delimiter','\t','precision',16) % number of elements +dlmwrite(fname,nelem,'delimiter','\t','precision',16) % number of elements +% dlmwrite(fname,nvert_e,'-append','delimiter','\t','precision',16) % number of vertices per element + +% % dlmwrite(fname,xc(:),'-append','delimiter','\t','precision',16) % coordinates +% dlmwrite(fname,xc([1,2],:).','-append','delimiter','\t','precision',16) % coordinates + +% subdomain info +tmp=int32(ia(4,:)); +% % dlmwrite(fname,tmp(:),'-append','delimiter','\t','precision',16) % connectivity in Matlab indexing +% dlmwrite(fname,tmp(:,:).','-append','delimiter','\t','precision',16) % connectivity in Matlab indexing +dlmwrite(fname,tmp(:,:).','-append','delimiter','\t') % connectivity in Matlab indexing + +% % store only start and end point of boundary edges, +% nbedges = size(e,2); +% dlmwrite(fname,nbedges,'-append','delimiter','\t','precision',16) % number boundary edges +% tmp=int32(e(1:2,:)); +% % dlmwrite(fname,tmp(:),'-append','delimiter','\t','precision',16) % boundary edges in Matlab indexing +% dlmwrite(fname,tmp(:,:).','-append','delimiter','\t','precision',16) % boundary edges in Matlab indexing + +end diff --git a/Sheet7/E14/jacob_template/cuthill_mckee_ordering.cpp b/Sheet7/E14/jacob_template/cuthill_mckee_ordering.cpp new file mode 100644 index 0000000..3aaad3f --- /dev/null +++ b/Sheet7/E14/jacob_template/cuthill_mckee_ordering.cpp @@ -0,0 +1,218 @@ +//======================================================================= +// Copyright 1997, 1998, 1999, 2000 University of Notre Dame. +// Authors: Andrew Lumsdaine, Lie-Quan Lee, Jeremy G. Siek +// +// This file is part of the Boost Graph Library +// +// You should have received a copy of the License Agreement for the +// Boost Graph Library along with the software; see the file LICENSE. +// If not, contact Office of Research, University of Notre Dame, Notre +// Dame, IN 46556. +// +// Permission to modify the code and to distribute modified code is +// granted, provided the text of this NOTICE is retained, a notice that +// the code was modified is included with the above COPYRIGHT NOTICE and +// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE +// file is distributed with the modified code. +// +// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. +// By way of example, but not limitation, Licensor MAKES NO +// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY +// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS +// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS +// OR OTHER RIGHTS. + +// Modyfied 2019 by Gundolf Haase, University of Graz +//======================================================================= +#include "cuthill_mckee_ordering.h" + +#include +#include +#include +#include +#include +#include +#include + +/* + Sample Output + original bandwidth: 8 + Reverse Cuthill-McKee ordering starting at: 6 + 8 3 0 9 2 5 1 4 7 6 + bandwidth: 4 + Reverse Cuthill-McKee ordering starting at: 0 + 9 1 4 6 7 2 8 5 3 0 + bandwidth: 4 + Reverse Cuthill-McKee ordering: + 0 8 5 7 3 6 4 2 1 9 + bandwidth: 4 + */ + /* +int main(int, char *[]) +{ + using namespace boost; + using namespace std; + typedef adjacency_list > > Graph; + typedef graph_traits::vertex_descriptor Vertex; + typedef graph_traits::vertices_size_type size_type; + + typedef std::pair Pair; + Pair edges[14] = { Pair(0, 3), //a-d + Pair(0, 5), //a-f + Pair(1, 2), //b-c + Pair(1, 4), //b-e + Pair(1, 6), //b-g + Pair(1, 9), //b-j + Pair(2, 3), //c-d + Pair(2, 4), //c-e + Pair(3, 5), //d-f + Pair(3, 8), //d-i + Pair(4, 6), //e-g + Pair(5, 6), //f-g + Pair(5, 7), //f-h + Pair(6, 7) + }; //g-h + + Graph G(10); + for (int i = 0; i < 14; ++i) + add_edge(edges[i].first, edges[i].second, G); + + graph_traits::vertex_iterator ui, ui_end; + + property_map::type deg = get(vertex_degree, G); + for (boost::tie(ui, ui_end) = vertices(G); ui != ui_end; ++ui) + deg[*ui] = degree(*ui, G); + + property_map::type + index_map = get(vertex_index, G); + + std::cout << "original bandwidth: " << bandwidth(G) << std::endl; + + std::vector inv_perm(num_vertices(G)); + std::vector perm(num_vertices(G)); + { + Vertex s = vertex(6, G); + //reverse cuthill_mckee_ordering + cuthill_mckee_ordering(G, s, inv_perm.rbegin(), get(vertex_color, G), + get(vertex_degree, G)); + cout << "Reverse Cuthill-McKee ordering starting at: " << s << endl; + cout << " "; + for (std::vector::const_iterator i = inv_perm.begin(); + i != inv_perm.end(); ++i) + cout << index_map[*i] << " "; + cout << endl; + + for (size_type c = 0; c != inv_perm.size(); ++c) + perm[index_map[inv_perm[c]]] = c; + std::cout << " bandwidth: " + << bandwidth(G, make_iterator_property_map(&perm[0], index_map, perm[0])) + << std::endl; + } + { + Vertex s = vertex(0, G); + //reverse cuthill_mckee_ordering + cuthill_mckee_ordering(G, s, inv_perm.rbegin(), get(vertex_color, G), + get(vertex_degree, G)); + cout << "Reverse Cuthill-McKee ordering starting at: " << s << endl; + cout << " "; + for (std::vector::const_iterator i = inv_perm.begin(); + i != inv_perm.end(); ++i) + cout << index_map[*i] << " "; + cout << endl; + + for (size_type c = 0; c != inv_perm.size(); ++c) + perm[index_map[inv_perm[c]]] = c; + std::cout << " bandwidth: " + << bandwidth(G, make_iterator_property_map(&perm[0], index_map, perm[0])) + << std::endl; + } + + { + //reverse cuthill_mckee_ordering + cuthill_mckee_ordering(G, inv_perm.rbegin(), get(vertex_color, G), + make_degree_map(G)); + + cout << "Reverse Cuthill-McKee ordering:" << endl; + cout << " "; + for (std::vector::const_iterator i = inv_perm.begin(); + i != inv_perm.end(); ++i) + cout << index_map[*i] << " "; + cout << endl; + + for (size_type c = 0; c != inv_perm.size(); ++c) + perm[index_map[inv_perm[c]]] = c; + std::cout << " bandwidth: " + << bandwidth(G, make_iterator_property_map(&perm[0], index_map, perm[0])) + << std::endl; + } + return 0; +} +*/ + + +// ------------- Modifications by Gundolf Haase +// std::vector _edges; //!< edges of mesh (vertices ordered ascending) + +using namespace boost; +using namespace std; +typedef adjacency_list > > Graph; +typedef graph_traits::vertex_descriptor Vertex; +typedef graph_traits::vertices_size_type size_type; + +typedef std::pair Pair; + +vector cuthill_mckee_reordering(vector const &_edges) +{ + size_t const nnodes = *max_element(cbegin(_edges), cend(_edges)) + 1; + cout << "NNODES = " << nnodes << endl; + //size_t const nedges = _edges.size()/2; + + Graph G(nnodes); + for (size_t i = 0; i < _edges.size(); i+=2) + add_edge(_edges[i], _edges[i+1], G); + + graph_traits::vertex_iterator ui, ui_end; + + property_map::type deg = get(vertex_degree, G); + for (boost::tie(ui, ui_end) = vertices(G); ui != ui_end; ++ui) + deg[*ui] = degree(*ui, G); + + property_map::type + index_map = get(vertex_index, G); + + std::cout << "original bandwidth: " << bandwidth(G) << std::endl; + + std::vector inv_perm(num_vertices(G)); + //std::vector perm(num_vertices(G)); + std::vector perm(num_vertices(G)); + + { + Vertex s = vertex(nnodes/2, G); + //reverse cuthill_mckee_ordering + cuthill_mckee_ordering(G, s, inv_perm.rbegin(), get(vertex_color, G), + get(vertex_degree, G)); + //cout << "Reverse Cuthill-McKee ordering starting at: " << s << endl; + //cout << " "; + //for (std::vector::const_iterator i = inv_perm.begin(); i != inv_perm.end(); ++i) + //cout << index_map[*i] << " "; + //cout << endl; + + for (size_type c = 0; c != inv_perm.size(); ++c) + perm[index_map[inv_perm[c]]] = c; + std::cout << "improved bandwidth: " + << bandwidth(G, make_iterator_property_map(&perm[0], index_map, perm[0])) + << std::endl; + } + + assert(perm.size()==nnodes); + + return perm; +} + +// ------------- end Modifications + + diff --git a/Sheet7/E14/jacob_template/cuthill_mckee_ordering.h b/Sheet7/E14/jacob_template/cuthill_mckee_ordering.h new file mode 100644 index 0000000..6891726 --- /dev/null +++ b/Sheet7/E14/jacob_template/cuthill_mckee_ordering.h @@ -0,0 +1,8 @@ +#ifndef CUTHILL_MCKEE_ORDERING +#define CUTHILL_MCKEE_ORDERING + +#include + +std::vector cuthill_mckee_reordering(std::vector const &_edges); + +#endif