MeshTopology.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "MeshTopology.h"

#include "bio/Cell.h"
#include "bio/Node.h"
#include "bio/Mesh.h"

#include <memory>
#include <string>
#include <vector>

using namespace std;
using namespace SimPT_Sim::Util;

namespace SimPT_Sim {

CSRMatrix MeshTopology::NodeCellNodeIncidence(shared_ptr<Mesh> mesh)
{
        const auto& nodes         = mesh->GetNodes();
        const size_t num_nodes    = nodes.size();

        // Nodes - cells incidence matrix
        CSRMatrix A;

        // The "running" A.row_ptr; I.e: start idx of a row (n_i) in A.col_ind
        size_t n_offset = 0;

        // For each node n_i, gather the cell indices (c_i) the node belongs to
        for (size_t n_i = 0; n_i < num_nodes; ++n_i) {
                const auto node = nodes[n_i];
                // Begin a new row by storing its start offset
                A.row_ptr.push_back(n_offset);
                for (const auto nbr : mesh->GetNodeOwningNeighbors(node)) {
                        const int c_i = nbr.GetCell()->GetIndex();
                        // I assume nodes sharing the boundary polygon are independent
                        if (c_i != -1) {
                                // Put a 1 at A[n_i, c_i] only if it wasn't there yet.
                                // (i.e: if c_i is NOT in [A.col_ind[A.row_ptr[n_i]], ... end])
                                if (find(A.col_ind.begin() + A.row_ptr.back(), A.col_ind.end(), c_i) == A.col_ind.end()) {
                                        A.col_ind.push_back(c_i);
                                        ++n_offset;
                                }
                        }
                }
        }

        // Add the last element to A.row_ptr, defined as length(col_ind)
        // See documentation of CSRMatrix.
        A.row_ptr.push_back(A.col_ind.size());
        return A;
}

CSRMatrix MeshTopology::NodeEdgeNodeIncidence(shared_ptr<Mesh> mesh)
{
        const auto& nodes        = mesh->GetNodes();
        const size_t num_nodes   = nodes.size();

        // Nodes - nodes incidence matrix
        CSRMatrix A;

        // The "running" A.row_ptr; I.e: start idx of a row (n_i) in A.col_ind
        size_t n_offset = 0;

        // For each node n_i, gather the indices of nodes it shares an edge with
        for (size_t n_i = 0; n_i < num_nodes; ++n_i) {
                const auto node = nodes[n_i];
                // Begin a new row by storing its start offset
                A.row_ptr.push_back(n_offset);
                for (const auto nbr : mesh->GetNodeOwningNeighbors(node)) {
                        const int n1_i = nbr.GetNb1()->GetIndex();
                        const int n2_i = nbr.GetNb2()->GetIndex();
                        // Put a 1 at A[n_i, n1/2_i] only if it wasn't there yet.
                        // (i.e: if n1/2_i is NOT in [A.col_ind[A.row_ptr[n_i]], ... end])
                        if (find(A.col_ind.begin() + A.row_ptr.back(), A.col_ind.end(), n1_i) == A.col_ind.end()) {
                                A.col_ind.push_back(n1_i);
                                ++n_offset;
                        }
                        if (find(A.col_ind.begin() + A.row_ptr.back(),
                                         A.col_ind.end(),
                                         n2_i) == A.col_ind.end()) {
                                A.col_ind.push_back(n2_i);
                                ++n_offset;
                        }
                }
        }

        // Add the last element to A.row_ptr, defined as length(col_ind)
        // cf. documentation of BinaryCSRMatrix.
        A.row_ptr.push_back(A.col_ind.size());
        return A;
}

} // namespace