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

#include "Cell.h"
#include "Edge.h"
#include "Node.h"
#include "NodeAttributes.h"
#include "WallAttributes.h"
#include "math/array3.h"
#include "math/constants.h"
#include "util/container/CircularIterator.h"
#include "util/container/ConstCircularIterator.h"

#include <functional>
#include <ostream>

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

namespace SimPT_Sim {

Wall::Wall(unsigned int index, Node* n1, Node* n2, Cell* c1, Cell* c2, unsigned int num_chem)
        : WallAttributes(num_chem),
          m_c1(c1), m_c2(c2), m_n1(n1), m_n2(n2), m_wall_index(index), m_length(0.0), m_length_dirty(true)
{
        assert( (n1 != nullptr) && "Wall::Wall> Attempting to build a node with nullptr cell!" );
        assert( (n2 != nullptr) && "Wall::Wall> Attempting to build a node with nullptr cell!" );
        assert( (c1 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr cell!" );
        assert( (c2 != nullptr) && "Wall::Wall> Attempting to build a wall with nullptr cell!" );
        assert( (c1 != c2) && "Wall::Wall> Attempting to build a wall between identical cells!" );

        std::array<double, 3> ref{{1., 0., 0.}};
        std::array<double, 3> tot1 = *n1;
        std::array<double, 3> tot2 = *n2;
        std::array<double, 3> tot3 = tot1 - tot2;
        double t2 = SignedAngle(tot3, ref);
        if ( ( t2 > (-pi() * 3 / 4) && t2 <= (-pi() / 4) ) || ( t2 > (pi() / 4) && t2 <= (pi() * 3 / 4) ) )
        {
                SetStrength(1.);
        } else {
                SetStrength(5.);//DDV: should be 1 for Blad, 5 for Wortel.
        }
}

void Wall::CalculateLength() const
{
        const auto& nodes    = m_c1->GetNodes();
        const auto  n_begin  = nodes.begin();
        const auto  n_end    = nodes.end();
        auto from = make_const_circular(nodes, std::find(n_begin, n_end, m_n1));
        auto to   = make_const_circular(nodes, std::find(n_begin, n_end, m_n2));

        assert(from != m_c1->GetNodes().end() && "N1 wall endpoint not in cell C1!");
        assert(to   != m_c1->GetNodes().end() && "N2 wall endpoint not in cell C1!");

        m_length = 0.0;
        do {
                m_length += Norm(**from - **std::next(from));
        } while (++from != to); // Loop until the next node is end point of wall
}

vector<Edge> Wall::GetEdges() const
{
        vector<Edge> v;
        const auto& nodes    = GetC1()->GetNodes();
        const auto cit_start = find(nodes.begin(), nodes.end(), GetN1());
        const auto cit_stop  = find(nodes.begin(), nodes.end(), GetN2());

        auto cit = make_const_circular(nodes, cit_start);
        do {
                v.push_back(Edge(*cit, *next(cit)));
        } while(++cit != cit_stop);
        return v;
}

double Wall::GetLength() const
{
        if (m_length_dirty) {
                CalculateLength();
                m_length_dirty = false;
        }
        return m_length;
}

array<double, 3> Wall::GetInfluxVector(Cell* c) const
{
        assert( (m_c1 == c || m_c2 == c) && "Wall::GetWallVector called with wrong cell!" );
        const double sign = ( c == m_c1 ? 1.0 : -1.0 );
        return sign * Orthogonalize(Normalize(*m_n2 - *m_n1));
}

array<double, 3> Wall::GetWallVector(Cell* c) const
{
        assert( (m_c1 == c || m_c2 ==c) && "Wall::GetWallVector called with wrong cell!" );
        const double sign = ( c == m_c1 ? 1.0 : -1.0 );
        return sign * (*m_n2 - *m_n1);
}

bool Wall::IsAtBoundary() const
{
        return (m_c1->GetIndex() == -1 || m_c2->GetIndex() == -1);
}

bool Wall::IsSam() const
{
        return GetN1()->IsSam() || GetN2()->IsSam();
}

bool Wall::IncludesEdge(const Edge& edge) const
{
        bool status = false;

        const auto& nodes     = GetC1()->GetNodes();
        const auto cit_start  = find(nodes.begin(), nodes.end(), GetN1());
        const auto cit_stop   = find(nodes.begin(), nodes.end(), GetN2());

        auto cit = make_const_circular(nodes, cit_start);
        do {
                status = (Edge(*cit, *next(cit)) == edge);
        } while(!status && ++cit != cit_stop);

        return status;
}

bool Wall::IncludesNode(const Node* node) const
{
        bool status = false;

        const auto& nodes     = GetC1()->GetNodes();
        const auto cit_start  = find(nodes.begin(), nodes.end(), GetN1());
        const auto cit_stop   = find(nodes.begin(), nodes.end(), GetN2());

        auto cit = make_const_circular(nodes, cit_start);
        do {
                status = (*cit == node);
        } while(!status && cit++ != cit_stop);

        return status;
}

bool Wall::IsIntersectingWithDivisionPlane(const array<double, 3>& p1, const array<double, 3>& p2) const
{
        const double x1 = (*m_n1)[0];
        const double y1 = (*m_n1)[1];
        const double x2 = (*m_n2)[0];
        const double y2 = (*m_n2)[1];
        const double x3 = p1[0];
        const double y3 = p1[1];
        const double x4 = p2[0];
        const double y4 = p2[1];
        const double ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3))
                        / ((y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1));

        // If ua is between 0 and 1, line p1 to p2 intersects the wall segment.
        return (ua >= 0. && ua <= 1.) ? true : false;

}

ostream& Wall::Print(ostream& os) const
{
        os << "Wall: " << GetIndex() << " "<< "{ " << m_n1->GetIndex() << "->" << m_n2->GetIndex()
                << ", " << m_c1->GetIndex() << " | " << m_c2->GetIndex() << "} ";
        return os;
}

void Wall::ReadPtree(boost::property_tree::ptree const& wall_pt)
{
        WallAttributes::ReadPtree(wall_pt.get_child("attributes"));
}

ptree Wall::ToPtree() const
{
        ptree ret;
        ret.put("id", GetIndex());
        ret.put("c1", m_c1->GetIndex());
        ret.put("c2", m_c2->GetIndex());
        ret.put("n1", m_n1->GetIndex());
        ret.put("n2", m_n2->GetIndex());
        ret.put("length", GetLength());
        ret.add_child("attributes", WallAttributes::ToPtree());

        return ret;
}

ostream& operator<<(ostream& os, const Wall& w)
{
        w.Print(os);
        return os;
}

}