TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
TRUST_2_MED.cpp
1/****************************************************************************
2* Copyright (c) 2026, 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 <TRUST_2_MED.h>
17#include <TRUSTTabs.h>
18#include <Char_ptr.h>
19#include <sys/stat.h>
20#include <Motcle.h>
21#include <fstream>
22#include <Noms.h>
23
24#include <medcoupling++.h>
25#ifdef MEDCOUPLING_
26#include <MCAuto.hxx>
27#include <MEDLoader.hxx>
28#endif
29
30#include <numeric>
31
32void traite_nom_fichier_med(Nom& nom_fic)
33{
34 Nom nom_fic2(nom_fic);
35 nom_fic2.prefix(".med");
36 if (nom_fic2==nom_fic)
37 {
38 Cerr<<"Error : the med field must have a .med extension"<<finl;
39 Cerr<<nom_fic2 <<" "<<nom_fic<<finl;
40 Process::exit();
41 }
42
43 nom_fic2.prefix("_0000");
44 nom_fic2+=".med";
45 nom_fic=nom_fic2.nom_me(nom_fic.me());
46 // on essaye en premier les fichiers _0000.med
47 {
48 std::ifstream test(nom_fic);
49 if (!test)
50 nom_fic=nom_fic2;
51 }
52 {
53 std::ifstream test(nom_fic);
54 if (!test)
55 {
56 // on essaye Cas_0000.med
57 nom_fic=nom_fic2.nom_me(0);
58 std::ifstream test2(nom_fic);
59 if (!test2)
60 {
61 Cerr<<"med file "<<nom_fic<<" not found."<<finl;
62 Process::exit();
63 }
64 }
65 }
66}
67
68extern "C" int MEDimport(char*,char*);
69void test_version(Nom& nom)
70{
71#ifdef MED_
72 // on regarde si le fichier est d'une version differente, si oui
73 // on cree un fichier au format MED majeur courant, et on change le nom du fichier
74 med_bool med_ok, hdf_ok;
75 if (MEDfileCompatibility(nom, &med_ok, &hdf_ok))
76 {
77 Cerr<<"Problem when trying to open the file "<<nom<<finl;
78 Process::exit();
79 }
80
81 if (hdf_ok && med_ok) return; //pas besoin de convertir
82
83 // On serialise pour eviter que le fichier soit cree plusieurs fois en //
84
85 Nom nom2bis("Conv_");
86 // ajout prefixe Conv_
87 Nom nom2 ;
88
89 char* nomtmp = new char[strlen(nom)+1];
90 int compteur=(int)strlen(nom);
91
92 strcpy(nomtmp,nom);
93
94 const char* ptr=nomtmp+compteur;
95 bool pas_de_slach=false;
96
97 while((*ptr!='/')&&(compteur>0)) // recherche du dernier slach dans le nom du fichier
98 {
99 ptr--;
100 compteur--;
101 if (*ptr=='/') pas_de_slach=true;
102 }
103
104 if (!pas_de_slach) // Cas ou le fichier med est dans le repertoire courant
105 {
106 nom2=nom2bis;
107 nom2+=nom;
108 }
109 else // Cas ou le fichier med est defini par un chemin
110 {
111 Nom cible(ptr);
112 cible.suffix("/"); // nom du fichier
113 nom2bis+= cible; // nom du fichier + prefix Conv_
114
115 Nom nomentier(nom);
116 nomentier.prefix(cible); // nom du chemin pour trouver fichier
117
118 nomentier+=nom2bis; // nom du chemin + prefixe Conv . nom fichier
119 nom2=nomentier;
120 }
121
122 delete [] nomtmp;
123
124 for (int p=0; p<Process::nproc(); p++)
125 {
126 if (p==Process::me())
127 {
128 bool a_creer= false;
129 std::ifstream test(nom2);
130 if (test)
131 {
132 struct stat org,newf;
133 stat(nom,&org);
134 stat(nom2,&newf);
135 if (org.st_mtime>newf.st_mtime)
136 a_creer=true;
137 }
138 else
139 a_creer=true;
140 if (a_creer)
141 {
142 Cerr<<"Creation of the file "<<nom2<<" to newer format."<<finl;
143 MEDimport(Char_ptr(nom),Char_ptr(nom2));
144 }
145 else
146 Cerr<<"The file "<<nom2<<" is up to date."<<finl;
147 }
148 Process::barrier();
149 }
150 nom=nom2;
151 return;
152#endif
153}
154
155void read_med_field_names(const Nom& nom_fic, Noms& noms_chps, ArrOfDouble& temps_sauv)
156{
157#if defined(MEDCOUPLING_) && defined(MED_)
158 using namespace std;
159 using namespace MEDCoupling;
160 using MCTimeLabel = pair< pair<int,int>, double>;
161 const string fnam(nom_fic.getString());
162
163 vector<string> nams(GetAllFieldNames(fnam));
164 noms_chps.resize((int)nams.size());
165 for(const auto& s: nams)
166 noms_chps.add(Nom(s));
167
168 vector<double> tims;
169 for(const auto& fldName: nams)
170 {
171 vector< MCTimeLabel > it(GetAllFieldIterations(fnam, fldName));
172 vector<double> other_tims;
173 other_tims.resize(it.size());
174 // return only the double representing the time:
175 auto lambd = [](const MCTimeLabel& cplx_pair) { return cplx_pair.second; };
176 transform(it.begin(), it.end(),
177 back_inserter(other_tims),
178 lambd);
179 std::sort(other_tims.begin(), other_tims.end());
180 if (tims.size() == 0)
181 tims = other_tims;
182 else if (tims != other_tims)
183 Process::exit("LireMED::read_med_field_names() - different timesteps between the various fields found in the MED file !!");
184 }
185 temps_sauv.resize((int)tims.size());
186 std::copy(tims.begin(), tims.end(), temps_sauv.addr());
187#else
188 med_non_installe();
189#endif
190}
191
192#ifdef MED_
193// renvoit le type med a partir du type trio
194med_geometry_type type_geo_trio_to_type_med(const Nom& type_elem_i,med_axis_type& rep)
195{
196 rep=MED_CARTESIAN;
197
198 // Strip any '_64_ prefix:
199 Motcle type_elem_0 = type_elem_i;
200 type_elem_0.prefix("_64");
201
202 // Check for axi:
203 Motcle type_elem = type_elem_0;
204 type_elem.prefix("_AXI");
205 if (type_elem != Motcle(type_elem_0))
206 {
207 rep=MED_SPHERICAL;
208 Cerr<<"#"<<type_elem<<"#"<<Motcle(type_elem_0)<<"#"<<(type_elem!=Motcle(type_elem_0))<<finl;
209 if (type_elem=="QUADRILATERE_2D")
210 type_elem="SEGMENT_2D";
211 if (type_elem=="RECTANGLE_2D")
212 {
213 type_elem="RECTANGLE";
214 rep=MED_CYLINDRICAL;
215 }
216 }
217 med_geometry_type type_elem_med;
218 if ((type_elem=="RECTANGLE") || (type_elem=="QUADRANGLE"))
219 type_elem_med=MED_QUAD4;
220 else if ((type_elem=="HEXAEDRE")|| (type_elem=="HEXAEDRE_VEF"))
221 type_elem_med=MED_HEXA8;
222 else if (type_elem=="TRIANGLE")
223 type_elem_med=MED_TRIA3;
224 else if (type_elem=="TETRAEDRE")
225 type_elem_med=MED_TETRA4;
226 else if ((type_elem=="SEGMENT_2D") || (type_elem=="SEGMENT"))
227 type_elem_med=MED_SEG2;
228 else if (type_elem=="TRIANGLE_3D")
229 type_elem_med=MED_TRIA3;
230 else if (type_elem=="QUADRANGLE_3D")
231 type_elem_med=MED_QUAD4;
232 else if (type_elem=="PRISME")
233 type_elem_med=MED_PENTA6;
234 else if (type_elem=="POLYEDRE")
235 type_elem_med=MED_POLYHEDRON;
236 else if (type_elem=="POLYGONE")
237 type_elem_med=MED_POLYGON;
238 else if (type_elem=="POLYGONE_3D")
239 type_elem_med=MED_POLYGON;
240 else if(type_elem=="SEGMENT")
241 type_elem_med=MED_SEG2;
242 else if(type_elem=="PRISME_HEXAG")
243 type_elem_med=MED_OCTA12;
244 else if(type_elem=="POINT_1D")
245 type_elem_med=MED_POINT1;
246 else if(type_elem=="POINT")
247 type_elem_med=MED_POINT1;
248 else
249 {
250 Cerr<<type_elem<< " no available code" <<finl;
251 Process::exit();
252 type_elem_med=MED_QUAD4;
253 }
254 return type_elem_med;
255}
256
257med_geometry_type type_geo_trio_to_type_med(const Nom& type_elem)
258{
259 med_axis_type rep;
260 return type_geo_trio_to_type_med(type_elem,rep);
261}
262
263#ifdef MEDCOUPLING_
264/*! @brief Return MEDCoupling type from TRUST type
265 * See file NormalizedGeometricTypes in MEDCoupling includes.
266 */
267INTERP_KERNEL::NormalizedCellType type_geo_trio_to_type_medcoupling(const Nom& type_elem_, int& mesh_dimension)
268{
269 Motcle type_elem;
270 type_elem=type_elem_;
271 type_elem.prefix("_AXI");
272 if (type_elem!=Motcle(type_elem_))
273 {
274 if (type_elem == "QUADRILATERE_2D")
275 type_elem = "SEGMENT_2D";
276 if (type_elem == "RECTANGLE_2D")
277 {
278 type_elem = "RECTANGLE";
279 }
280 }
281 mesh_dimension = -1;
282 INTERP_KERNEL::NormalizedCellType type_cell;
283 if ((type_elem=="RECTANGLE") || (type_elem=="QUADRANGLE") || (type_elem=="QUADRANGLE_3D"))
284 {
285 type_cell = INTERP_KERNEL::NORM_QUAD4;
286 mesh_dimension = 2;
287 }
288 else if ((type_elem=="HEXAEDRE") || (type_elem=="HEXAEDRE_VEF"))
289 {
290 type_cell = INTERP_KERNEL::NORM_HEXA8;
291 mesh_dimension = 3;
292 }
293 else if ((type_elem=="TRIANGLE") || (type_elem=="TRIANGLE_3D"))
294 {
295 type_cell = INTERP_KERNEL::NORM_TRI3;
296 mesh_dimension = 2;
297 }
298 else if (type_elem=="TETRAEDRE")
299 {
300 type_cell = INTERP_KERNEL::NORM_TETRA4;
301 mesh_dimension = 3;
302 }
303 else if ((type_elem=="SEGMENT") || (type_elem=="SEGMENT_2D"))
304 {
305 type_cell = INTERP_KERNEL::NORM_SEG2;
306 mesh_dimension = 1;
307 }
308 else if (type_elem=="PRISME")
309 {
310 type_cell = INTERP_KERNEL::NORM_PENTA6;
311 mesh_dimension = 3;
312 }
313 else if (type_elem=="POLYEDRE")
314 {
315 type_cell = INTERP_KERNEL::NORM_POLYHED;
316 mesh_dimension = 3;
317 }
318 else if ((type_elem=="POLYGONE") || (type_elem=="POLYGONE_3D"))
319 {
320 type_cell = INTERP_KERNEL::NORM_POLYGON;
321 mesh_dimension = 2;
322 }
323 else if(type_elem=="PRISME_HEXAG")
324 {
325 type_cell = INTERP_KERNEL::NORM_HEXGP12;
326 mesh_dimension = 3;
327 }
328 else if ((type_elem=="POINT") || (type_elem=="POINT_1D"))
329 {
330 type_cell = INTERP_KERNEL::NORM_POINT1;
331 mesh_dimension = 0;
332 }
333 else
334 {
335 Cerr<<type_elem<< " no available cell." <<finl;
336 Process::exit();
337 return INTERP_KERNEL::NORM_POINT1;
338 }
339 assert(mesh_dimension>=0);
340 return type_cell;
341}
342
343void fill_connectivity_from_mc_mesh(const MEDCoupling::MEDCouplingUMesh * mc_mesh, IntTab& face_sommet, IntTab& elem_face)
344{
345 using namespace MEDCoupling;
346 using DAI = MCAuto<DataArrayIdType>;
347
348 DAI desc(DataArrayIdType::New()), descIndx(DataArrayIdType::New()), revDesc(DataArrayIdType::New()), revDescIndx(DataArrayIdType::New());
349 MCAuto<MEDCoupling::MEDCouplingUMesh> mc_desc(mc_mesh->buildDescendingConnectivity(desc, descIndx, revDesc, revDescIndx));
350
351 DataArrayIdType* c(mc_desc->getNodalConnectivity()), *cI(mc_desc->getNodalConnectivityIndex());
352 const mcIdType *cP(c->getConstPointer()), *cIP(cI->getConstPointer());
353
354 // Fill face_sommet
355 int nb_faces = Process::check_int_overflow(mc_desc->getNumberOfCells());
356 DAI dsi = cI->deltaShiftIndex();
357 const mcIdType* dsiP = dsi->getConstPointer();
358 mcIdType idx_dnu;
359 int max_pts = static_cast<int>(dsi->getMaxValue(idx_dnu));
360 face_sommet.resize(nb_faces, max_pts-1);
361 face_sommet = -1;
362
363 for (auto i = 0; i < nb_faces; i++)
364 {
365 auto nb_pts = dsiP[i] - 1;
366 for (auto j = 0; j < nb_pts; j++)
367 face_sommet(i,j) = static_cast<int>(cP[cIP[i]+1+j]); // +1 to skip type information (NORM_xxx)
368 }
369
370 // Filling elem_face
371 int nb_elem = Process::check_int_overflow(mc_mesh->getNumberOfCells());
372 DAI dsi2 = descIndx->deltaShiftIndex();
373 int max_nb_of_fac = static_cast<int>(dsi2->getMaxValue(idx_dnu));
374 elem_face.resize(nb_elem, max_nb_of_fac);
375 elem_face = -1;
376 const mcIdType *descP(desc->getConstPointer()), *descIndxP(descIndx->getConstPointer());
377
378 for (auto i = 0; i < nb_elem; i++)
379 for (auto j = 0; j < descIndxP[i+1]-descIndxP[i]; j++)
380 elem_face(i, j) = static_cast<int>(descP[descIndxP[i]+j]);
381}
382
383#endif /* MEDCOUPLING_ */
384#endif /* MED_ */
385
386/*! @brief Passage de la connectivite TRUST a MED si toMED=true de MED a trio si toMED=false
387 */
388template <typename _SIZE_>
389void conn_trust_to_med(IntTab_T<_SIZE_>& les_elems, const Nom& type_elem, bool toMED)
390{
391#ifdef MED_
392 using int_t = _SIZE_;
393 using IntTab_t = IntTab_T<_SIZE_>;
394
395 int_t nele=les_elems.dimension(0);
396 // cas face_bord vide
397 if (nele==0) return;
398 med_geometry_type type_elem_med;
399 type_elem_med=type_geo_trio_to_type_med(type_elem);
400 IntTab_t les_elemsn(les_elems);
401 ArrOfInt filter;
402 switch (type_elem_med)
403 {
404 case MED_POINT1 :
405 filter.resize_array(1) ;
406 filter[0] = 0 ;
407 break ;
408 case MED_SEG2 :
409 filter.resize_array(2) ;
410 std::iota(filter.addr(), filter.addr()+2, 0);
411 break ;
412 case MED_SEG3 :
413 break ;
414 case MED_TRIA3 :
415 filter.resize_array(3) ;
416 std::iota(filter.addr(), filter.addr()+3, 0);
417 break ;
418 case MED_QUAD4 :
419 filter.resize_array(4) ;
420 filter[0] = 0 ;
421 filter[1] = 2 ;
422 filter[2] = 3 ;
423 filter[3] = 1 ;
424 break ;
425 case MED_TRIA6 :
426 case MED_QUAD8 :
427 break ;
428 case MED_TETRA4 :
429 filter.resize_array(4) ;
430 std::iota(filter.addr(), filter.addr()+4, 0);
431 break ;
432 case MED_PYRA5 :
433 filter.resize_array(5) ;
434 filter[0] = 0 ;
435 filter[1] = 3 ; // 2nd element in med are 4th in vtk (array begin at 0 !)
436 filter[2] = 2 ;
437 filter[3] = 1 ; // 4th element in med are 2nd in vtk (array begin at 0 !)
438 filter[4] = 4 ;
439 break ;
440 case MED_PENTA6 :
441 filter.resize_array(6) ;
442 std::iota(filter.addr(), filter.addr()+6, 0);
443 break ;
444 case MED_HEXA8 :
445 filter.resize_array(8) ;
446 filter[0] = 0 ;
447 filter[1] = 2 ;
448 filter[2] = 3 ;
449 filter[3] = 1 ;
450 filter[4] = 4 ;
451 filter[5] = 6 ;
452 filter[6] = 7 ;
453 filter[7] = 5 ;
454 break ;
455 case MED_TETRA10 :
456 case MED_PYRA13 :
457 case MED_PENTA15 :
458 case MED_HEXA20 :
459 break ;
460
461 case MED_OCTA12:
462 case MED_POLYGON:
463 case MED_POLYHEDRON:
464 {
465 int nb_som_max=les_elems.dimension_int(1);
466 filter.resize_array(nb_som_max);
467 std::iota(filter.addr(), filter.addr()+nb_som_max, 0);
468 break ;
469 }
470 default:
471 Cerr<<"case not scheduled"<<finl;
472 Process::exit();
473 }
474
475 int ns=filter.size_array();
476 if (ns==0)
477 {
478 Cerr<<"Problem for filtering operation "<<finl;
479 Process::exit();
480 }
481 if (toMED)
482 {
483 for (int_t el=0; el<nele; el++)
484 for (int n=0; n<ns; n++)
485 les_elems(el,n)=les_elemsn(el,filter[n]);
486 }
487 else
488 {
489 for (int_t el=0; el<nele; el++)
490 for (int n=0; n<ns; n++)
491 les_elems(el,filter[n])=les_elemsn(el,n);
492 }
493#endif
494}
495
496// Explicit instanciation
497template void conn_trust_to_med(IntTab_T<int>& les_elems, const Nom& type_elem, bool toMED);
498#if INT_is_64_ == 2
499template void conn_trust_to_med(IntTab_T<trustIdType>& les_elems, const Nom& type_elem, bool toMED);
500#endif
Char_ptr
class Char_ptr Une chaine de caractere pour nommer les objets de TRUST
Definition Char_ptr.h:28
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Nom::nom_me
Nom nom_me(int, const char *prefix=0, int without_padding=0) const
Insere _prefix000n (n=me() ou nproc()) dans un nom de fichier (par ex:toto.
Definition Nom.cpp:387
Nom::prefix
Nom & prefix(const char *const)
Definition Nom.cpp:329
Nom::getString
const std::string & getString() const
Definition Nom.h:92
Noms
Un tableau de chaine de caracteres (VECT(Nom)).
Definition Noms.h:26
Process::check_int_overflow
static int check_int_overflow(trustIdType)
Definition Process.cpp:428
Process::nproc
static int nproc()
renvoie le nombre de processeurs dans le groupe courant Voir Comm_Group::nproc() et PE_Groups::curren...
Definition Process.cpp:104
Process::barrier
static void barrier()
Synchronise tous les processeurs du groupe courant (attend que tous les processeurs soient arrives a ...
Definition Process.cpp:136
Process::me
static int me()
renvoie mon rang dans le groupe de communication courant.
Definition Process.cpp:125
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
TRUSTArray::size_array
_SIZE_ size_array() const
Definition TRUSTArray.tpp:187
TRUSTArray::addr
_TYPE_ * addr()
Definition TRUSTArray.tpp:159
TRUSTArray::resize_array
void resize_array(_SIZE_ new_size, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTArray.tpp:43
TRUSTArray::resize
void resize(_SIZE_ new_size, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTArray.h:156
TRUSTTab::dimension_int
int dimension_int(int d) const
Definition TRUSTTab.tpp:152
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133