TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
OrienteFacesBord.cpp
1/****************************************************************************
2* Copyright (c) 2024, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <OrienteFacesBord.h>
17#include <Motcle.h>
18
19#include <SFichier.h>
20
21Implemente_instanciable_32_64(OrienteFacesBord_32_64,"OrienteFacesBord",Interprete_geometrique_base_32_64<_T_>);
22// XD orientefacesbord interprete orientefacesbord INHERITS_BRACE Keyword to modify the order of the boundary vertices
23// XD_CONT included in a domain, such that the surface normals are outer pointing.
24// XD attr domain_name ref_domaine domain_name REQ Name of domain.
25
26
27template <typename _SIZE_>
28Sortie& OrienteFacesBord_32_64<_SIZE_>::printOn(Sortie& os) const
29{
30 return Interprete::printOn(os);
31}
32
33template <typename _SIZE_>
34Entree& OrienteFacesBord_32_64<_SIZE_>::readOn(Entree& is)
35{
36 return Interprete::readOn(is);
37}
38
39template <typename _SIZE_>
40Entree& OrienteFacesBord_32_64<_SIZE_>::interpreter_(Entree& is)
41{
42 this->associer_domaine(is);
43 oriente_faces_bord(this->domaine().le_nom());
44 return is;
45}
46
47
48template <typename _SIZE_>
49void OrienteFacesBord_32_64<_SIZE_>::oriente_faces_bord(const Nom& nom_dom)
50{
51 Domaine_t& dom=ref_cast(Domaine_t, this->objet(nom_dom));
52 const Nom type=dom.type_elem()->que_suis_je();
53 constexpr double e=0.001;
54 SmallArrOfTID_t elem(1);
55 DoubleTab pos(1,this->dimension);
56 const DoubleTab_t& coords=dom.coord_sommets();
57 const int nombre_faces_rayonnantes=dom.nb_bords()+dom.nb_raccords();
58 for (int iface=0; iface<nombre_faces_rayonnantes; iface++)
59 {
60 const IntTab_t& som=(iface<dom.nb_bords()?dom.bord(iface).faces().les_sommets():dom.raccord(iface-dom.nb_bords())->faces().les_sommets());
61 const int_t nb_faces=som.dimension(0);
62 IntTab_t new_faces(nb_faces,som.dimension_int(1));
63
64 for (int_t j=0; j<nb_faces; j++)
65 {
66 if (this->dimension==2)
67 {
68 if (this->axi)
69 {
70 const double r0=coords(som(j,0),0);
71 const double t0=coords(som(j,0),1);
72 const double r1=coords(som(j,1),0);
73 double t1=coords(som(j,1),1);
74 if (t1<t0) t1+=2*M_PI;
75 if (est_egal(r0,r1))
76 {
77 pos(0,1)=0.5*(t0+t1);
78 if (t0<t1)
79 pos(0,0)=(1-e)*r0;
80 else
81 pos(0,0)=(1+e)*r0;
82 }
83 else if (est_egal(t0,t1))
84 {
85 pos(0,0)=0.5*(r0+r1);
86 if (r0<r1)
87 pos(0,1)=(1+e)*t0;
88 else
89 pos(0,1)=(1-e)*t0;
90 }
91 }
92 else
93 {
94 const double x0=coords(som(j,0),0);
95 const double y0=coords(som(j,0),1);
96 const double x1=coords(som(j,1),0);
97 const double y1=coords(som(j,1),1);
98 const double nx=y0-y1;
99 const double ny=x1-x0;
100 const double alpha=10*Objet_U::precision_geom/(sqrt(nx*nx+ny*ny));
101 pos(0,0)=0.5*(x0+x1)+alpha*nx;
102 pos(0,1)=0.5*(y0+y1)+alpha*ny;
103 }
104 dom.chercher_elements(pos,elem);
105 if (elem(0)!=-1)
106 {
107 new_faces(j,0) = som(j,0);
108 new_faces(j,1) = som(j,1);
109 }
110 else
111 {
112 new_faces(j,0) = som(j,1);
113 new_faces(j,1) = som(j,0);
114 }
115 }
116 else
117 {
118 if (this->axi)
119 {
120 const double r0=coords(som(j,0),0);
121 const double t0=coords(som(j,0),1);
122 const double z0=coords(som(j,0),2);
123 const double r3=coords(som(j,3),0);
124 const double z3=coords(som(j,3),2);
125 double t3=coords(som(j,3),1);
126 // Cas ou on repasse a teta=0;
127 if (t3<t0) t3+=2*M_PI;
128
129 if (est_egal(r0,r3))
130 {
131 pos(0,1)=0.5*(t0+t3);
132 if (t0<t3)
133 pos(0,0)=(1+e)*r0;
134 else
135 pos(0,0)=(1-e)*r0;
136 pos(0,2)=0.5*(z0+z3);
137 }
138 else if (est_egal(z0,z3))
139 {
140 pos(0,0)=0.5*(r0+r3);
141 pos(0,1)=0.5*(t0+t3);
142 if (r0<r3)
143 pos(0,2)=(1+e)*z0;
144 else
145 pos(0,2)=(1-e)*z0;
146 }
147 else if (est_egal(t0,t3))
148 {
149 pos(0,0)=0.5*(r0+r3);
150 pos(0,2)=0.5*(z0+z3);
151 if (r0<r3)
152 pos(0,1)=(1+e)*t0;
153 else
154 pos(0,1)=(1-e)*t0;
155 }
156 dom.chercher_elements(pos,elem);
157 if (elem(0)!=-1)
158 {
159 new_faces(j,0) = som(j,0);
160 new_faces(j,1) = som(j,2);
161 new_faces(j,2) = som(j,3);
162 new_faces(j,3) = som(j,1);
163 }
164 else
165 {
166 new_faces(j,0) = som(j,0);
167 new_faces(j,1) = som(j,1);
168 new_faces(j,2) = som(j,3);
169 new_faces(j,3) = som(j,2);
170 }
171 }
172 else
173 {
174 const double x0=coords(som(j,0),0);
175 const double y0=coords(som(j,0),1);
176 const double z0=coords(som(j,0),2);
177 const double x1=coords(som(j,1),0);
178 const double y1=coords(som(j,1),1);
179 const double z1=coords(som(j,1),2);
180 const double x2=coords(som(j,2),0);
181 const double y2=coords(som(j,2),1);
182 const double z2=coords(som(j,2),2);
183 const double nx=(y1-y0)*(z2-z0)-(z1-z0)*(y2-y0);
184 const double ny=(z1-z0)*(x2-x0)-(x1-x0)*(z2-z0);
185 const double nz=(x1-x0)*(y2-y0)-(y1-y0)*(x2-x0);
186 const double alpha=100*Objet_U::precision_geom/(sqrt(nx*nx+ny*ny+nz*nz));
187 pos(0,0)=(x0+x1+x2)/3.+alpha*nx;
188 pos(0,1)=(y0+y1+y2)/3.+alpha*ny;
189 pos(0,2)=(z0+z1+z2)/3.+alpha*nz;
190 dom.chercher_elements(pos,elem);
191 /*
192 if (x0>0 && y0>0 && x1>0 && y1>0 && x2>0 && y2>0)
193 {
194 Cout << "nx = " << nx << " " << ny << " " << nz << finl;
195 Cout << "pos = " << pos << finl;
196 Cout << elem(0) << finl;
197 } */
198 if (elem(0)!=-1)
199 {
200 //Cout << "non inversion" << finl;
201 if (type=="Tetraedre")
202 {
203 new_faces(j,0) = som(j,0);
204 new_faces(j,1) = som(j,1);
205 new_faces(j,2) = som(j,2);
206
207 }
208 else
209 {
210 new_faces(j,0) = som(j,0);
211 new_faces(j,1) = som(j,1);
212 new_faces(j,2) = som(j,3);
213 new_faces(j,3) = som(j,2);
214 }
215 }
216 else
217 {
218 //Cout << "inversion" << finl;
219 if (type=="Tetraedre")
220 {
221 new_faces(j,0) = som(j,0);
222 new_faces(j,1) = som(j,2);
223 new_faces(j,2) = som(j,1);
224 }
225 else
226 {
227 new_faces(j,0) = som(j,0);
228 new_faces(j,1) = som(j,2);
229 new_faces(j,2) = som(j,3);
230 new_faces(j,3) = som(j,1);
231 }
232 }
233 if (elem(0)!=-1)
234 {
235 const int_t el=elem(0);
236 pos(0,0)=(x0+x1+x2)/3.-alpha*nx;
237 pos(0,1)=(y0+y1+y2)/3.-alpha*ny;
238 pos(0,2)=(z0+z1+z2)/3.-alpha*nz;
239 dom.chercher_elements(pos,elem);
240 if (elem(0)!=-1)
241 {
242 Cerr << "Case not implemented in the OrienteFacesBord algorithm! " << finl;
243 Cerr << "The face " << j << " is surrounded by the two elements " << el << " " << elem(0) << finl;
244 Cerr << alpha*nx << finl;
245 Cerr << alpha*ny << finl;
246 Cerr << alpha*nz << finl;
247 Cerr << "Contact TRUST support." <<finl;
248 this->exit();
249 }
250 }
251 }
252 }
253 }
254 if (iface<dom.nb_bords())
255 dom.bord(iface).faces().les_sommets().ref(new_faces);
256 else
257 dom.raccord(iface-dom.nb_bords())->faces().les_sommets().ref(new_faces);
258 }
259
260}
261
262template class OrienteFacesBord_32_64<int>;
263#if INT_is_64_ == 2
264template class OrienteFacesBord_32_64<trustIdType>;
265#endif
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Interprete_geometrique_base_32_64
classe Interprete_geometrique_base .
Definition Interprete_geometrique_base.h:37
Interprete_geometrique_base_32_64::associer_domaine
void associer_domaine(Nom &nom_dom)
Definition Interprete_geometrique_base.cpp:32
Interprete_geometrique_base_32_64::domaine
Domaine_t & domaine(int i=0)
Definition Interprete_geometrique_base.h:48
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Objet_U::Entree
friend class Entree
Definition Objet_U.h:76
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::precision_geom
static double precision_geom
Definition Objet_U.h:86
Objet_U::le_nom
virtual const Nom & le_nom() const
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
Definition Objet_U.cpp:319
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
OrienteFacesBord_32_64
class OrienteFacesBord
Definition OrienteFacesBord.h:33
OrienteFacesBord_32_64::interpreter_
Entree & interpreter_(Entree &) override
Definition OrienteFacesBord.cpp:40
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52