NodeInserter.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 "NodeInserter.h"

#include "bio/Cell.h"
#include "bio/Edge.h"
#include "bio/Node.h"
#include "bio/Mesh.h"
#include "bio/Wall.h"
#include "sim/Sim.h"
#include "util/misc/Exception.h"
#include "util/misc/log_debug.h"

using namespace std;
using namespace boost::property_tree;
using namespace SimPT_Sim;
using namespace SimPT_Sim::Util;

namespace SimPT_Sim{

NodeInserter::NodeInserter()
        : m_mesh(nullptr), m_target_node_distance(0.0), m_yielding_threshold(0.0)
{
}

NodeInserter::NodeInserter(const CoreData& cd)
        : m_mesh(nullptr), m_target_node_distance(0.0), m_yielding_threshold(0.0)
{
        Initialize(cd);
}


void NodeInserter::Initialize(const CoreData& cd)
{
        try {
                assert( cd.Check() && "CoreData not ok.");
                m_cd     = cd;
                m_mesh   = m_cd.m_mesh.get();

                m_target_node_distance  = m_cd.m_parameters->get<double>("cell_mechanics.target_node_distance");
                m_yielding_threshold    = m_cd.m_parameters->get<double>("cell_mechanics.yielding_threshold");
        }
        catch(exception& e) {
                const string here = string(VL_HERE) + " exception:\n";
                throw Exception(here + e.what());
        }
}

void NodeInserter::InsertNode(Edge& edge)
{
        Wall *wall = m_mesh->FindWall(edge);

        // Construct a new node in the middle of the edge
        const array<double, 3> pos = (*edge.GetFirst() + *edge.GetSecond()) / 2.0;
        auto new_node = m_mesh->BuildNode(pos, m_mesh->IsAtBoundary(&edge));
        new_node->SetFixed(edge.GetFirst()->IsFixed() && edge.GetSecond()->IsFixed());
        new_node->SetSam(new_node->IsAtBoundary() &&
                (edge.GetFirst()->IsSam() || edge.GetSecond()->IsSam()));

        UpdateNeighbors(edge, new_node);
        m_mesh->GetNodeOwningWalls(new_node) = { wall };
}

unsigned int NodeInserter::InsertNodes()
{
        unsigned int insert_count = 0U;
        const auto& walls = m_mesh->GetWalls();

        for (const auto& wall : walls) {
                if ( wall->GetC1()->IsFixed() || wall->GetC2()->IsFixed()) {
                        continue;
                }
                vector<Edge> edges = wall->GetEdges();
                for (auto edge : edges) {
                        const double l = Norm(*edge.GetFirst() - *edge.GetSecond());
                        if (l > m_yielding_threshold * m_target_node_distance) {
                                InsertNode(edge);
                                insert_count++;
                        }
                }
        }

        return insert_count;
}

void NodeInserter::UpdateNeighbors(Edge& edge, Node* new_node)
{
        list<NeighborNodes> owners;

        // Lambdas  to compare neighbors.
        auto cmp = [](const NeighborNodes& n1, const NeighborNodes& n2) {
                return n1.GetCell()->GetIndex() < n2.GetCell()->GetIndex();
        };
        auto eql = [](const NeighborNodes& n1, const NeighborNodes& n2) {
                return n1.GetCell()->GetIndex() == n2.GetCell()->GetIndex();
        };

        // Push all cells owning the two nodes of the divided edges onto a list.
        auto& owners_first = m_mesh->GetNodeOwningNeighbors(edge.GetFirst());
        copy(owners_first.begin(), owners_first.end(), back_inserter(owners));
        auto& owners_second = m_mesh->GetNodeOwningNeighbors(edge.GetSecond());
        copy(owners_second.begin(), owners_second.end(), back_inserter(owners));
        owners.sort(cmp);

        // The duplicates in this list indicate cells owning this edge.
        auto c = owners.begin();
        while (c != owners.end()) {
                c = adjacent_find(c, owners.end(), eql);

                if (c != owners.end()) {
                        m_mesh->AddNodeToCell(c->GetCell(), new_node, edge.GetFirst(), edge.GetSecond());
                } else {
                        break;
                }
                c++;
        }

        // Repair neighborhood lists in a second loop, to make sure it is up to date
        while (c != owners.end()) {
                c = adjacent_find(c, owners.end(), eql);
                // repair neighborhood lists of cell and Wall lists
                if (!c->GetCell()->IsBoundaryPolygon()) {
                        m_mesh->ConstructNeighborList(c->GetCell());
                }
                c++;
        }
}

} // namespace