| #ifndef _generate_matrix_structure_hpp_ |
| #define _generate_matrix_structure_hpp_ |
| |
| //@HEADER |
| // ************************************************************************ |
| // |
| // MiniFE: Simple Finite Element Assembly and Solve |
| // Copyright (2006-2013) Sandia Corporation |
| // |
| // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive |
| // license for use of this work by or on behalf of the U.S. Government. |
| // |
| // This library is free software; you can redistribute it and/or modify |
| // it under the terms of the GNU Lesser General Public License as |
| // published by the Free Software Foundation; either version 2.1 of the |
| // License, or (at your option) any later version. |
| // |
| // This library is distributed in the hope that it will be useful, but |
| // WITHOUT ANY WARRANTY; without even the implied warranty of |
| // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| // Lesser General Public License for more details. |
| // |
| // You should have received a copy of the GNU Lesser General Public |
| // License along with this library; if not, write to the Free Software |
| // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 |
| // USA |
| // |
| // ************************************************************************ |
| //@HEADER |
| |
| #include <sstream> |
| #include <stdexcept> |
| #include <map> |
| #include <algorithm> |
| |
| #include <simple_mesh_description.hpp> |
| #include <SparseMatrix_functions.hpp> |
| #include <box_utils.hpp> |
| #include <utils.hpp> |
| |
| #ifdef HAVE_MPI |
| #include <mpi.h> |
| #endif |
| |
| namespace miniFE { |
| |
| template<typename MatrixType> |
| int |
| generate_matrix_structure(const simple_mesh_description<typename MatrixType::GlobalOrdinalType>& mesh, |
| MatrixType& A) |
| { |
| int myproc = 0; |
| #ifdef HAVE_MPI |
| MPI_Comm_rank(MPI_COMM_WORLD, &myproc); |
| #endif |
| |
| int threw_exc = 0; |
| try { |
| |
| typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; |
| typedef typename MatrixType::LocalOrdinalType LocalOrdinal; |
| |
| int global_nodes_x = mesh.global_box[0][1]+1; |
| int global_nodes_y = mesh.global_box[1][1]+1; |
| int global_nodes_z = mesh.global_box[2][1]+1; |
| Box box; |
| copy_box(mesh.local_box, box); |
| |
| //num-owned-nodes in each dimension is num-elems+1 |
| //only if num-elems > 0 in that dimension *and* |
| //we are at the high end of the global range in that dimension: |
| if (box[0][1] > box[0][0] && box[0][1] == mesh.global_box[0][1]) ++box[0][1]; |
| if (box[1][1] > box[1][0] && box[1][1] == mesh.global_box[1][1]) ++box[1][1]; |
| if (box[2][1] > box[2][0] && box[2][1] == mesh.global_box[2][1]) ++box[2][1]; |
| |
| GlobalOrdinal global_nrows = global_nodes_x; |
| global_nrows *= global_nodes_y*global_nodes_z; |
| |
| GlobalOrdinal nrows = get_num_ids<GlobalOrdinal>(box); |
| try { |
| A.reserve_space(nrows, 27); |
| } |
| catch(std::exception& exc) { |
| std::ostringstream osstr; |
| osstr << "One of A.rows.resize, A.row_offsets.resize, A.packed_cols.reserve or A.packed_coefs.reserve: nrows=" <<nrows<<": "; |
| osstr << exc.what(); |
| std::string str1 = osstr.str(); |
| throw std::runtime_error(str1); |
| } |
| |
| std::vector<GlobalOrdinal> rows(nrows); |
| std::vector<LocalOrdinal> row_offsets(nrows+1); |
| std::vector<int> row_coords(nrows*3); |
| |
| unsigned roffset = 0; |
| GlobalOrdinal nnz = 0; |
| |
| for(int iz=box[2][0]; iz<box[2][1]; ++iz) { |
| for(int iy=box[1][0]; iy<box[1][1]; ++iy) { |
| for(int ix=box[0][0]; ix<box[0][1]; ++ix) { |
| GlobalOrdinal row_id = |
| get_id<GlobalOrdinal>(global_nodes_x, global_nodes_y, global_nodes_z, |
| ix, iy, iz); |
| rows[roffset] = mesh.map_id_to_row(row_id); |
| row_coords[roffset*3] = ix; |
| row_coords[roffset*3+1] = iy; |
| row_coords[roffset*3+2] = iz; |
| row_offsets[roffset++] = nnz; |
| |
| for(int sz=-1; sz<=1; ++sz) { |
| for(int sy=-1; sy<=1; ++sy) { |
| for(int sx=-1; sx<=1; ++sx) { |
| GlobalOrdinal col_id = |
| get_id<GlobalOrdinal>(global_nodes_x, global_nodes_y, global_nodes_z, |
| ix+sx, iy+sy, iz+sz); |
| if (col_id >= 0 && col_id < global_nrows) { |
| ++nnz; |
| } |
| } |
| } |
| } |
| |
| } |
| } |
| } |
| |
| if (roffset != nrows) { |
| throw std::runtime_error("ERROR in generate_matrix_structure, roffset != nrows."); |
| } |
| row_offsets[roffset] = nnz; |
| |
| init_matrix(A, rows, row_offsets, row_coords, |
| global_nodes_x, global_nodes_y, global_nodes_z, global_nrows, mesh); |
| } |
| catch(...) { |
| std::cout << "proc " << myproc << " threw an exception in generate_matrix_structure, probably due to running out of memory." << std::endl; |
| threw_exc = 1; |
| } |
| #ifdef HAVE_MPI |
| int global_throw = 0; |
| MPI_Allreduce(&threw_exc, &global_throw, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
| threw_exc = global_throw; |
| #endif |
| if (threw_exc) { |
| return 1; |
| } |
| |
| return 0; |
| } |
| |
| }//namespace miniFE |
| |
| #endif |
| |