Code_TYMPAN  4.4.0
Industrial site acoustic simulation
PostTreatment.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) <2012> <EDF-R&D> <FRANCE>
3  * This program is free software; you can redistribute it and/or modify
4  * it under the terms of the GNU General Public License as published by
5  * the Free Software Foundation; either version 2 of the License, or
6  * (at your option) any later version.
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
10  * See the GNU General Public License for more details.
11  * You should have received a copy of the GNU General Public License along
12  * with this program; if not, write to the Free Software Foundation, Inc.,
13  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
14  */
15 
16 #include <set>
17 #include <vector>
18 #include <map>
19 
21 #include "Geometry/Triangle.h"
22 #include "Geometry/Cylindre.h"
23 #include "Tools/Logger.h"
24 #include "PostTreatment.h"
25 
26 #include <math.h>
27 
28 typedef std::pair<unsigned int, unsigned int> segment;
29 typedef std::map<segment, std::vector<Shape*>> mapSegmentShapes;
30 
32 {
33 
34  Triangle* triangle = dynamic_cast<Triangle*>(shape);
35  if (!shape)
36  {
37  return;
38  }
39 
40  std::vector<unsigned int>* vertices = triangle->getLocalVertices();
41  for (unsigned int i = 0; i < vertices->size(); i++)
42  {
43  segment seg1 = std::pair<unsigned int, unsigned int>(vertices->at(i),
44  vertices->at((i + 1) % (vertices->size())));
45  mapSegmentShapes::iterator it = currentMap.find(seg1);
46 
47  if (it == currentMap.end())
48  {
49  std::vector<Shape*> newList;
50  newList.push_back(shape);
51  currentMap.insert(std::pair<segment, std::vector<Shape*>>(seg1, newList));
52  }
53  else
54  {
55  it->second.push_back(shape);
56  }
57 
58  segment seg2 = std::pair<unsigned int, unsigned int>(vertices->at((i + 1) % (vertices->size())),
59  vertices->at(i));
60  it = currentMap.find(seg2);
61 
62  if (it == currentMap.end())
63  {
64  std::vector<Shape*> newList;
65  newList.push_back(shape);
66  currentMap.insert(std::pair<segment, std::vector<Shape*>>(seg2, newList));
67  }
68  else
69  {
70  it->second.push_back(shape);
71  }
72  }
73 }
74 
75 bool isAcceptableEdge(const segment& seg, Shape* p1, Shape* p2, decimal& angleOuverture,
76  vector<bool>& isGround)
77 {
78 
79  // First we test is the current segment does not belong to both a construct and the ground
80  if (!p1->isSol() && !p2->isSol() && isGround[seg.first] && isGround[seg.second])
81  {
82  return false;
83  }
84 
85  /*
86  we test if p1 vs p2 angle is greater than PI/2+angleMax
87  else they are considered as colinear and no diffraction cylinder is built
88  seg
89  |
90  p1 |p2/
91  | /
92  __________|/
93  comp \
94  \ normal
95  \
96 
97  */
98 
99  // Minimal angle (other PI) between two face to allow building of a diffraction cylinder
100  float angleMax = (decimal)(AcousticRaytracerConfiguration::get()->AngleDiffMin * M_PI / 180);
101 
102  // Compute "mean" normal between the two faces
103  vec3 normal = p1->getNormal() + p2->getNormal();
104  normal.normalize();
105 
106  // Search for a vertex of the first shape not belonging to seg
107  std::vector<unsigned int>* vertices = p1->getLocalVertices();
108  vec3 otherVertex;
109  for (unsigned int i = 0; i < vertices->size(); i++)
110  {
111  if (vertices->at(i) != seg.first && vertices->at(i) != seg.second)
112  {
113  otherVertex = p1->getVertices()->at(vertices->at(i));
114  break;
115  }
116  }
117 
118  // Find a vector othogonal to seg and lying in the first shape
119  vec3 proj =
120  otherVertex.closestPointOnLine(p1->getVertices()->at(seg.first), p1->getVertices()->at(seg.second));
121  vec3 comp = otherVertex - proj;
122  comp.normalize();
123 
124  // Compute the angle betwwen comp and the "mean" normal
125  decimal cos_angle = comp.dot(normal);
126  angleOuverture = 2.f * acos(cos_angle);
127 
128  if (cos_angle < cos(M_PI / 2. + angleMax / 2.))
129  {
130  return true;
131  }
132 
133  return false;
134 }
135 
137 {
138  // define diffraction cylinder diameter
139  float cylinderThick = (float)AcousticRaytracerConfiguration::get()->CylindreThick;
140 
141  // Create a list of segments common to two faces
142  mapSegmentShapes segmentList;
143  std::vector<Shape*>* shapes = scene->getShapes();
144  std::vector<bool> isGround(scene->getVertices()->size(), false);
145 
146  for (unsigned int i = 0; i < shapes->size(); i++)
147  {
148  Shape* shape = shapes->at(i);
149  registerSegmentsFromShapes(shape, segmentList);
150 
151  if (shape->isSol())
152  {
153  for (unsigned int j = 0; j < shape->getLocalVertices()->size(); j++)
154  {
155  isGround[shape->getLocalVertices()->at(j)] = true;
156  }
157  }
158  }
159 
160  ss << segmentList.size() << " segments ont ete enregistres." << std::endl;
161 
162  std::set<segment> validSegment;
163 
164  for (mapSegmentShapes::iterator it = segmentList.begin(); it != segmentList.end(); it++)
165  {
166  // If segment is owned by 2 shapes, we test angle between them
167  if (it->second.size() == 2)
168  {
169  std::set<segment>::iterator itset1 = validSegment.find(it->first);
170  std::set<segment>::iterator itset2 =
171  validSegment.find(std::pair<unsigned int, unsigned int>(it->first.second, it->first.first));
172  decimal angleOuverture = NAN;
173  if (itset1 == validSegment.end() && itset2 == validSegment.end() &&
174  isAcceptableEdge(it->first, it->second.at(0), it->second.at(1), angleOuverture, isGround))
175  {
176  validSegment.insert(it->first);
177  Cylindre* cylindre = new Cylindre(it->second.at(0), it->second.at(1), scene->getVertices(),
178  it->first.first, it->first.second, cylinderThick);
179  cylindre->setAngleOuverture(angleOuverture);
180  scene->addShape(cylindre);
181  ss << "Ajout du segment (" << it->first.first << "," << it->first.second << ")" << std::endl;
182  }
183  }
184  }
185 
186  ss << validSegment.size() << " cylindre ont ete ajoutes." << std::endl;
187 
188  return true;
189 }
190 
191 #ifdef _ALLOW_TARGETING_
192 bool PostTreatment::findTargetsForNMPB(Scene* scene, std::vector<Recepteur>& recepteurs,
193  TargetManager& targetManager, decimal density)
194 {
195 
196  std::vector<Shape*>* shapes = scene->getShapes();
197  for (unsigned int i = 0; i < shapes->size(); i++)
198  {
199  if (dynamic_cast<Triangle*>(shapes->at(i)))
200  {
201  /*Triangle* t = dynamic_cast<Triangle*>(shapes->at(i));
202  if(!(t->getMaterial()->isNatural)){
203  vec3 normale = t->getNormal();
204  vec2 normale2D = vec2(sqrt(normale.x * normale.x + normale.y * normale.y),normale.z);
205  normale2D.normalize();
206  //Si l'angle de la normale par rapport a l'horizontale est plus faible que 15 degres pour etre
207  coherent avec la definition de surface verticale de la NMPB if(acos(normale2D.dot(vec2(1.,0.))) <
208  M_PI / 12.){ std::vector<vec3> samples; t->sample(density,samples);
209  targetManager.registerTargets(samples);
210  }
211  }*/
212  }
213  else if (dynamic_cast<Cylindre*>(shapes->at(i)))
214  {
215  Cylindre* c = dynamic_cast<Cylindre*>(shapes->at(i));
216  std::vector<vec3> samples;
217  c->sample(density, samples);
218  targetManager.registerTargets(samples);
219  }
220  }
221 
222  for (unsigned int i = 0; i < recepteurs.size(); i++)
223  {
224  targetManager.registerTarget(recepteurs.at(i).getPosition());
225  }
226 
227  return true;
228 }
229 
230 void PostTreatment::appendDirectionToSources(TargetManager* targets, std::vector<Source>& sources)
231 {
232  /*for(unsigned int i = 0; i < sources.size(); i++){
233  vec3 pos = sources.at(i).getPosition();
234  for(unsigned int j = 0; j < targets.size(); j++){
235  vec3 newDir = targets.at(j) - pos;
236  newDir.normalize();
237  sources.at(i).addDirection(newDir);
238  }
239  sources.at(i).setInitialRayCount(targets.size());
240  }*/
241  for (unsigned int i = 0; i < sources.size(); i++)
242  {
243  sources.at(i).setTargetManager(targets);
244  sources.at(i).setInitialTargetCount(sources.at(i).getInitialRayCount());
245  }
246 }
247 #endif
std::stringstream ss
Definition: Logger.cpp:21
NxReal c
Definition: NxVec3.cpp:317
std::pair< unsigned int, unsigned int > segment
bool isAcceptableEdge(const segment &seg, Shape *p1, Shape *p2, decimal &angleOuverture, vector< bool > &isGround)
void registerSegmentsFromShapes(Shape *shape, mapSegmentShapes &currentMap)
std::map< segment, std::vector< Shape * > > mapSegmentShapes
double CylindreThick
Diffraction cylinder diameter.
static AcousticRaytracerConfiguration * get()
Get access to the configuration.
Cylinder class.
Definition: Cylindre.h:27
void setAngleOuverture(decimal angle)
Definition: Cylindre.h:43
This class mainly define a mesh (list of Shape) used by the Simulation object.
Definition: Scene.h:51
std::vector< Shape * > * getShapes()
Return all the shapes.
Definition: Scene.h:88
void addShape(Shape *shape)
Add a shape to the list.
Definition: Scene.h:68
std::vector< vec3 > * getVertices()
Return all the vertices.
Definition: Scene.h:97
base class for shapes (Cylindre, Mesh, Sphere, Triangle,...)
Definition: Shape.h:57
vector< vec3 > * getVertices()
Definition: Shape.h:127
vector< unsigned int > * getLocalVertices()
Get local vertices.
Definition: Shape.h:133
bool isSol() const
Get/Set the flag _isSol (ground or not)
Definition: Shape.h:194
virtual vec3 getNormal(const vec3 pos=vec3())
Get normal.
Definition: Shape.h:145
Class to manage targets.
Definition: TargetManager.h:27
bool registerTargets(std::vector< vec3 > &newTargets)
Register a vector of targets.
bool registerTarget(const vec3 newTarget)
Register a new target.
Triangle class.
Definition: Triangle.h:25
#define M_PI
Pi.
Definition: color.cpp:25
bool constructEdge(Scene *scene)
Build the edges list of a scene.
float decimal
Definition: mathlib.h:45
base_vec3< decimal > vec3
Definition: mathlib.h:381