TYAcousticModel.cpp

Go to the documentation of this file.

/*
 * Copyright (C) <2012-2024> <EDF-DTG> <FRANCE>
 * This file is part of Code_TYMPAN (R).
 * Code_TYMPAN (R) is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * Code_TYMPAN (R) 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 General Public License for more details.
 * You should have received a copy of the GNU General Public License along
 * with Code_TYMPAN (R). If not, see <https://www.gnu.org/licenses/>.
 */
 
#include <deque>
#include <list>
#include <cmath>
#include <algorithm>
#include "Tympan/core/defines.h"
#include "Tympan/models/solver/config.h"
#include "Tympan/solvers/9613/TYTrajet.h"
#include "Tympan/solvers/9613/TYSolver.h"
#include "Tympan/models/common/mathlib.h"
#include "Tympan/geometric_methods/AcousticRaytracer/Geometry/Scene.h"
#include "TYAcousticModel.h"
 
#include "Tympan/solvers/9613/TYChemin.h"
 
TYAcousticModel::TYAcousticModel()
    : _useSol(true), _useReflex(false), _propaCond(0), _interference(false), _paramH(10.0)
{
}
 
TYAcousticModel::~TYAcousticModel() {}
 
void TYAcousticModel::init()
{
    tympan::LPSolverConfiguration config = tympan::SolverConfiguration::get();
    // Calcul avec sol reel
    _useSol = config->UseRealGround;
    // Calcul avec reflexion sur les parois verticales
    _useReflex = config->UseReflection;
    // Calcul en conditions favorables
    _propaCond = config->PropaConditions;
    // Definit l'atmosphere courante du site
    double pression = config->AtmosPressure;
    double temperature = config->AtmosTemperature;
    double hygrometrie = config->AtmosHygrometry;
 
    _pSolverAtmos =
        std::unique_ptr<AtmosphericConditions>(new AtmosphericConditions(pression, temperature, hygrometrie));
    // Calcul avec interference
    _interference = config->ModSummation;
    // Compute wave length
    computeWaveLength();
    // Coefficient multiplicateur pour le calcul des reflexions supplementaires en condition favorable
    _paramH = config->H1parameter;
}
 
OPlan TYAcousticModel::buildMeanSlopePlan(const OSegment3D& penteMoyenne) const
{
    OPoint3D pt1{penteMoyenne._ptA};
    OPoint3D pt2{penteMoyenne._ptB};
    OPoint3D pt3{};
    pt3._z = pt1._z;
 
    if (pt2._x != pt1._x)
    {
        pt3._y = pt1._y + 1;
        pt3._x = (pt1._y - pt2._y) * (pt3._y - pt1._y) / (pt2._x - pt1._x) + (pt1._x);
    }
    else
    {
        if (pt1._y != pt2._y)
        {
            pt3._x = pt1._x + 1;
            pt3._y = (pt2._x - pt1._x) * (pt3._x - pt1._x) / (pt1._y - pt2._y) + (pt1._y);
        }
        else // pt1 and pt2 coincide, we shift pt2, building an horizontal plane
        {
            pt2._x = pt1._x + 1;
            pt2._y = pt1._y - 1;
            pt3._y = pt1._y + 1;
            pt3._x = (pt1._y - pt2._y) * (pt3._y - pt1._y) / (pt2._x - pt1._x) + (pt1._x);
        }
    }
 
    return OPlan{pt1, pt2, pt3};
}
 
void TYAcousticModel::meanSlope(const OSegment3D& director, OSegment3D& slope) const
{
    // Search for primitives under the two segment extremities
 
    // To begin : initialize slope
    slope = director;
    TYSolver& solver = getSolver();
 
    // first one
    OPoint3D pt = director._ptA;
    pt._z += 1000.;
    vec3 start = OPoint3Dtovec3(pt);
    Ray ray1(start, vec3(0., 0., -1.));
 
    std::list<Intersection> LI;
 
    double distance1 = static_cast<double>(solver.getScene()->getAccelerator()->traverse(&ray1, LI));
    assert(distance1 > 0.);
    assert(!LI.empty());
 
    unsigned int indexFace = LI.begin()->p->getPrimitiveId();
 
    // Avoid cases where the extremities are above infrastructure elements
    while (solver.getTabPolygon()[indexFace].is_infra())
    {
        start.z = (decimal)std::min(std::min(solver.getTabPolygon()[indexFace].tabPoint[0]._z,
                                             solver.getTabPolygon()[indexFace].tabPoint[1]._z),
                                    solver.getTabPolygon()[indexFace].tabPoint[2]._z);
        Ray ray(start, vec3(0, 0, -1));
        ray.setMaxt(20000);
        std::list<Intersection> LI2;
        double distance = static_cast<double>(solver.getScene()->getAccelerator()->traverse(&ray, LI2));
        // assert(distance > 0.);
        if (LI2.empty() || distance < 0.)
            break;
        distance1 += distance;
        indexFace = LI2.begin()->p->getPrimitiveId();
    }
 
    // Second one
    LI.clear();
    pt = director._ptB;
    pt._z += 1000.;
    start = OPoint3Dtovec3(pt);
    Ray ray2(start, vec3(0., 0., -1.));
 
    double distance2 = static_cast<double>(solver.getScene()->getAccelerator()->traverse(&ray2, LI));
    // An error can occur if some elements are outside of the grip (emprise)
    assert(distance2 > 0.);
    assert(!LI.empty());
 
    indexFace = LI.begin()->p->getPrimitiveId();
 
    // Avoid cases where the extremities are above infrastructure elements
    while (solver.getTabPolygon()[indexFace].is_infra())
    {
        start.z = (decimal)std::min(std::min(solver.getTabPolygon()[indexFace].tabPoint[0]._z,
                                             solver.getTabPolygon()[indexFace].tabPoint[1]._z),
                                    solver.getTabPolygon()[indexFace].tabPoint[2]._z);
        Ray ray(start, vec3(0, 0, -1));
        ray.setMaxt(20000);
        std::list<Intersection> LI2;
        double distance = static_cast<double>(solver.getScene()->getAccelerator()->traverse(&ray, LI2));
        // assert(distance > 0.);
        if (LI2.empty() || distance < 0.)
            break;
        distance2 += distance;
        indexFace = LI2.begin()->p->getPrimitiveId();
    }
 
    // Compute projection on the ground of segment points suppose sol is under the points ...
    slope._ptA._z = director._ptA._z - (distance1 - 1000.);
    slope._ptB._z = director._ptB._z - (distance2 - 1000.);
}