TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
distances_EF.cpp
1/****************************************************************************
2* Copyright (c) 2023, 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#include <distances_EF.h>
16#include <Domaine.h>
17#include <Motcle.h>
18
19double norm_3D_vit1(const DoubleTab& vit,int fac,int num1,
20 const Domaine_EF& domaine,
21 double& val1,double& val2,double& val3)
22{
23 const DoubleTab& face_normale = domaine.face_normales();
24 // fac numero de la face a paroi fixe
25 double r0,r1,r2;
26 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
27 int nsom = domaine.nb_som_face();
28 ArrOfDouble vit_face(3);
29 vit_face=0.;
30 for(int comp=0; comp<3; comp++)
31 {
32 for(int jsom=0; jsom<nsom; jsom++)
33 {
34 int num_som = domaine.face_sommets(num1, jsom);
35 vit_face[comp]+=vit(num_som,comp);
36 }
37 }
38 double v1=vit_face[0]/nsom;
39 double v2=vit_face[1]/nsom;
40 double v3=vit_face[2]/nsom;
41 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
42
43 double psc = r0*v1+r1*v2+r2*v3;
44 // val1,val2 val3 sont les vitesses tangentielles
45 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
46 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
47 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
48
49 return norm_vit;
50}
51
52double norm_2D_vit1_k(const DoubleTab& vit,int fac,int num1,
53 const Domaine_EF& domaine,
54 double& val1,double& val2)
55{
56 const DoubleTab& face_normale = domaine.face_normales();
57 // fac numero de la face a paroi fixe
58 double r0,r1;
59 calcule_r0r1(face_normale,fac,r0,r1);
60 int nsom = domaine.nb_som_face();
61 ArrOfDouble vit_face(2);
62 vit_face=0.;
63 for(int comp=0; comp<2; comp++)
64 {
65 for(int jsom=0; jsom<nsom; jsom++)
66 {
67 int num_som = domaine.face_sommets(num1, jsom);
68 vit_face[comp]+=vit(num_som,comp);
69 }
70 }
71 double v1=vit_face[0]/nsom;
72 double v2=vit_face[1]/nsom;
73 double norm_vit = vitesse_tangentielle(v1,v2,r0,r1);
74
75 double psc = r0*v1+r1*v2;
76
77 // val1 et val2 sont les vitesses tangetentielles
78 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
79 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
80
81 return norm_vit;
82}
83
84double norm_3D_vit2(const DoubleTab& vit,int fac,int num1,
85 const Domaine_EF& domaine,
86 double& val1,double& val2,double& val3)
87{
88 const DoubleTab& face_normale = domaine.face_normales();
89 // fac numero de la face a paroi defilante
90 double r0,r1,r2;
91 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
92
93 int nsom = domaine.nb_som_face();
94 ArrOfDouble vit_face(3);
95 ArrOfDouble vit_face2(3);
96 vit_face=0.;
97 vit_face2=0.;
98 for(int comp=0; comp<3; comp++)
99 {
100 for(int jsom=0; jsom<nsom; jsom++)
101 {
102 int num_som = domaine.face_sommets(num1, jsom);
103 vit_face[comp]+=vit(num_som,comp);
104 int num_som2 = domaine.face_sommets(fac, jsom);
105 vit_face2[comp]+=vit(num_som2,comp);
106 }
107 }
108 double v1=(vit_face[0]-vit_face2[0])/nsom;
109 double v2=(vit_face[1]-vit_face2[1])/nsom;
110 double v3=(vit_face[2]-vit_face2[2])/nsom;
111 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
112 // val1,val2 val3 sont les vitesses tangentielles
113 double psc = r0*v1+r1*v2+r2*v3;
114 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
115 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
116 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
117
118 return norm_vit;
119}
120
121double norm_2D_vit2_k(const DoubleTab& vit,int fac,int num1,
122 const Domaine_EF& domaine,
123 double& val1,double& val2)
124{
125 const DoubleTab& face_normale = domaine.face_normales();
126 // fac numero de la face a paroi defilante
127 double r0,r1;
128 calcule_r0r1(face_normale,fac,r0,r1);
129
130 int nsom = domaine.nb_som_face();
131 ArrOfDouble vit_face(2);
132 ArrOfDouble vit_face2(2);
133 vit_face=0.;
134 vit_face2=0.;
135 for(int comp=0; comp<2; comp++)
136 {
137 for(int jsom=0; jsom<nsom; jsom++)
138 {
139 int num_som = domaine.face_sommets(num1, jsom);
140 vit_face[comp]+=vit(num_som,comp);
141 int num_som2 = domaine.face_sommets(fac, jsom);
142 vit_face2[comp]+=vit(num_som2,comp);
143 }
144 }
145 double v1=(vit_face[0]-vit_face2[0])/nsom;
146 double v2=(vit_face[1]-vit_face2[1])/nsom;
147 double norm_vit = vitesse_tangentielle(v1,v2,r0,r1);
148 // val1 et val2 sont les vitesses tangetentielles
149 double psc = r0*v1+r1*v2;
150 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
151 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
152
153 return norm_vit;
154}
155
156double distance_2D(int fac,int elem,const Domaine_EF& domaine)
157{
158 const DoubleTab& xp = domaine.xp(); // centre de gravite des elements
159 const DoubleTab& xv = domaine.xv(); // centre de gravite des faces
160
161 const DoubleTab& face_normale = domaine.face_normales();
162 double r0,r1;
163 calcule_r0r1(face_normale,fac,r0,r1);
164
165 double x0=xv(fac,0);
166 double y0=xv(fac,1);
167 double x1=xp(elem,0);
168 double y1=xp(elem,1);
169
170 return fabs(r0*(x1-x0)+r1*(y1-y0));
171}
172
173double distance_3D(int fac,int elem,const Domaine_EF& domaine)
174{
175 const DoubleTab& xp = domaine.xp(); // centre de gravite des elements
176 const DoubleTab& xv = domaine.xv(); // centre de gravite des faces
177 const DoubleTab& face_normale = domaine.face_normales();
178 double r0,r1,r2;
179 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
180 double x0=xv(fac,0);
181 double y0=xv(fac,1);
182 double z0=xv(fac,2);
183 double x1=xp(elem,0);
184 double y1=xp(elem,1);
185 double z1=xp(elem,2);
186
187 return fabs(r0*(x1-x0)+r1*(y1-y0)+r2*(z1-z0));
188}
189
190double norm_vit_lp(const ArrOfDouble& vit,int face,const Domaine_EF& domaine,ArrOfDouble& val)
191{
192
193 const DoubleTab& face_normale = domaine.face_normales();
194 int dim=Objet_U::dimension;
195 ArrOfDouble r(dim);
196 double psc,norm_vit;
197
198 for(int i=0; i<dim; i++) r[i]=face_normale(face,i);
199
200 r/=norme_array(r);
201 psc = dotproduct_array(r,vit);
202
203 if(dim==3) norm_vit = vitesse_tangentielle(vit[0],vit[1],vit[2],r[0],r[1],r[2]);
204 else norm_vit = vitesse_tangentielle(vit[0],vit[1],r[0],r[1]);
205
206 for(int i=0; i<dim; i++) val[i]=(vit[i]-psc*r[i])/(norm_vit+DMINFLOAT);
207
208 return norm_vit;
209}
210
211double norm_vit_lp_k(const DoubleTab& vit,int face,int face_b,const Domaine_EF& domaine,ArrOfDouble& val,int is_defilante)
212{
213 int dimension=Objet_U::dimension;
214 //assert(face_b<domaine.premiere_face_int()); assert a ameliorer ou faces_virt...
215 if (is_defilante==0)
216 if (dimension==3)
217 return norm_3D_vit1(vit,face_b,face,domaine,val[0],val[1],val[2]);
218 else
219 return norm_2D_vit1_k(vit,face_b,face,domaine,val[0],val[1]);
220
221 else
222 {
223 if (dimension==3)
224 return norm_3D_vit2(vit,face_b,face,domaine,val[0],val[1],val[2]);
225 else
226 return norm_2D_vit2_k(vit,face_b,face,domaine,val[0],val[1]);
227 }
228}
229
230double distance_face_elem(int fac,int elem,const Domaine_EF& domaine)
231{
232 int dimension=Objet_U::dimension;
233 if (dimension==3)
234 return distance_3D(fac,elem,domaine);
235 else
236 return distance_2D(fac,elem,domaine);
237}
238
Domaine_EF
class Domaine_EF
Definition Domaine_EF.h:59
Objet_U::dimension
static int dimension
Definition Objet_U.h:99