66/* */
77/*----------------------------------------------------------------------------*/
88/* */
9- /* Description: open a mesh sequentialy and read/write */
10- /* its fields in multithread */
9+ /* Description: open a mesh sequentialy then read & write its fields */
10+ /* in a multithread way. Also create and write an */
11+ /* associated solb in parallel */
1112/* Author: Loic MARECHAL */
1213/* Creation date: feb 25 2025 */
13- /* Last modification: mar 11 2026 */
14+ /* Last modification: jun 19 2026 */
1415/* */
1516/*----------------------------------------------------------------------------*/
1617
3132
3233typedef struct
3334{
34- int ver , dim , TetTyp , VerTyp , NmbCpu , (* TetTab )[5 ], * RefTab ;
35- int64_t NmbVer , NmbTet , InpMsh , OutMsh , ParIdx ;
36- char * InpNam , * OutNam ;
35+ int ver , dim , TetTyp , TriTyp , VerTyp , NmbCpu ;
36+ int (* TriTab )[4 ], (* TetTab )[5 ], * RefTab ;
37+ int64_t NmbVer , NmbTri , NmbTet , InpMsh , OutMsh , OutSol , ParIdx ;
38+ char * InpNam , * OutNam , * SolNam ;
3739 double (* VerTab )[3 ];
3840}MshSct ;
3941
4042
4143/*----------------------------------------------------------------------------*/
42- /* Read the vertices */
44+ /* Read the vertices in parallel */
4345/*----------------------------------------------------------------------------*/
4446
4547void ScaVer (int BegIdx , int EndIdx , int PthIdx , MshSct * msh )
@@ -53,7 +55,7 @@ void ScaVer(int BegIdx, int EndIdx, int PthIdx, MshSct *msh)
5355 exit (1 );
5456 }
5557
56- printf ("THREAD %3d: read vertices %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
58+ printf ("THREAD %3d: read vertices %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
5759 GmfGetBlock (InpMsh , GmfVertices , BegIdx , EndIdx , 0 , NULL , NULL ,
5860 GmfDoubleVec , 3 , msh -> VerTab [ BegIdx ], msh -> VerTab [ EndIdx ],
5961 GmfInt , & msh -> RefTab [ BegIdx ], & msh -> RefTab [ EndIdx ]);
@@ -63,7 +65,30 @@ void ScaVer(int BegIdx, int EndIdx, int PthIdx, MshSct *msh)
6365
6466
6567/*----------------------------------------------------------------------------*/
66- /* Read the tets */
68+ /* Read the triangles in parallel */
69+ /*----------------------------------------------------------------------------*/
70+
71+ void ScaTri (int BegIdx , int EndIdx , int PthIdx , MshSct * msh )
72+ {
73+ int ver , dim ;
74+ int64_t InpMsh ;
75+
76+ if (!(InpMsh = GmfOpenMesh (msh -> InpNam , GmfRead , & ver , & dim )))
77+ {
78+ printf ("Thread %d failed to reopen file %s\n" , PthIdx , msh -> InpNam );
79+ exit (1 );
80+ }
81+
82+ printf ("THREAD %3d: read triangles %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
83+ GmfGetBlock (InpMsh , GmfTriangles , BegIdx , EndIdx , 0 , NULL , NULL ,
84+ GmfIntVec , 4 , msh -> TriTab [ BegIdx ], msh -> TriTab [ EndIdx ]);
85+
86+ GmfCloseMesh (InpMsh );
87+ }
88+
89+
90+ /*----------------------------------------------------------------------------*/
91+ /* Read the tets in parallel */
6792/*----------------------------------------------------------------------------*/
6893
6994void ScaTet (int BegIdx , int EndIdx , int PthIdx , MshSct * msh )
@@ -77,7 +102,7 @@ void ScaTet(int BegIdx, int EndIdx, int PthIdx, MshSct *msh)
77102 exit (1 );
78103 }
79104
80- printf ("THREAD %3d: read tets %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
105+ printf ("THREAD %3d: read tets %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
81106 GmfGetBlock (InpMsh , GmfTetrahedra , BegIdx , EndIdx , 0 , NULL , NULL ,
82107 GmfIntVec , 5 , msh -> TetTab [ BegIdx ], msh -> TetTab [ EndIdx ]);
83108
@@ -86,7 +111,7 @@ void ScaTet(int BegIdx, int EndIdx, int PthIdx, MshSct *msh)
86111
87112
88113/*----------------------------------------------------------------------------*/
89- /* Write the vertices */
114+ /* Write the vertices in parallel */
90115/*----------------------------------------------------------------------------*/
91116
92117void RecVer (int BegIdx , int EndIdx , int PthIdx , MshSct * msh )
@@ -100,7 +125,7 @@ void RecVer(int BegIdx, int EndIdx, int PthIdx, MshSct *msh)
100125 exit (1 );
101126 }
102127
103- printf ("THREAD %3d: write vertices %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
128+ printf ("THREAD %3d: write vertices %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
104129 res = GmfSetBlock (OutMsh , GmfVertices , BegIdx , EndIdx , 0 , NULL , NULL ,
105130 GmfDoubleVec , 3 , msh -> VerTab [ BegIdx ], msh -> VerTab [ EndIdx ],
106131 GmfInt , & msh -> RefTab [ BegIdx ], & msh -> RefTab [ EndIdx ]);
@@ -110,7 +135,30 @@ void RecVer(int BegIdx, int EndIdx, int PthIdx, MshSct *msh)
110135
111136
112137/*----------------------------------------------------------------------------*/
113- /* Write the tets */
138+ /* Write the triangles in parallel */
139+ /*----------------------------------------------------------------------------*/
140+
141+ void RecTri (int BegIdx , int EndIdx , int PthIdx , MshSct * msh )
142+ {
143+ int res ;
144+ int64_t OutMsh ;
145+
146+ if (!(OutMsh = GmfOpenMesh (msh -> OutNam , GmfStartParallelWrite , msh -> ver , msh -> dim )))
147+ {
148+ printf ("Thread %d failed to reopen file %s\n" , PthIdx , msh -> OutNam );
149+ exit (1 );
150+ }
151+
152+ printf ("THREAD %3d: write triangles %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
153+ res = GmfSetBlock (OutMsh , GmfTriangles , BegIdx , EndIdx , 0 , NULL , NULL ,
154+ GmfIntVec , 4 , msh -> TriTab [ BegIdx ], msh -> TriTab [ EndIdx ]);
155+
156+ GmfCloseUnfinishedMesh (OutMsh );
157+ }
158+
159+
160+ /*----------------------------------------------------------------------------*/
161+ /* Write the tets in parallel */
114162/*----------------------------------------------------------------------------*/
115163
116164void RecTet (int BegIdx , int EndIdx , int PthIdx , MshSct * msh )
@@ -124,35 +172,67 @@ void RecTet(int BegIdx, int EndIdx, int PthIdx, MshSct *msh)
124172 exit (1 );
125173 }
126174
127- printf ("THREAD %3d: write tets %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
175+ printf ("THREAD %3d: write tets %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
128176 res = GmfSetBlock (OutMsh , GmfTetrahedra , BegIdx , EndIdx , 0 , NULL , NULL ,
129177 GmfIntVec , 5 , msh -> TetTab [ BegIdx ], msh -> TetTab [ EndIdx ]);
130178
131179 GmfCloseUnfinishedMesh (OutMsh );
132180}
133181
134182
183+ /*----------------------------------------------------------------------------*/
184+ /* Write the solutions at vertices in parallel */
185+ /*----------------------------------------------------------------------------*/
186+
187+ void RecSol (int BegIdx , int EndIdx , int PthIdx , MshSct * msh )
188+ {
189+ int res ;
190+ int64_t OutSol ;
191+ int TypTab [3 ], SizTab [3 ];
192+ char * BegTab [3 ], * EndTab [3 ];
193+
194+ TypTab [0 ] = GmfDoubleVec ;
195+ SizTab [0 ] = 3 ;
196+ BegTab [0 ] = (char * )msh -> VerTab [ BegIdx ];
197+ EndTab [0 ] = (char * )msh -> VerTab [ EndIdx ];
198+
199+ if (!(OutSol = GmfOpenMesh (msh -> SolNam , GmfStartParallelWrite , msh -> ver , msh -> dim )))
200+ {
201+ printf ("Thread %d failed to reopen file %s\n" , PthIdx , msh -> OutNam );
202+ exit (1 );
203+ }
204+
205+ printf ("THREAD %3d: write solution %10d -> %10d\n" , PthIdx , BegIdx , EndIdx );
206+
207+ res = GmfSetBlock (OutSol , GmfSolAtVertices , BegIdx , EndIdx , 0 , NULL , NULL ,
208+ GmfDoubleVec , 3 , msh -> VerTab [ BegIdx ], msh -> VerTab [ EndIdx ]);
209+
210+ GmfCloseUnfinishedMesh (OutSol );
211+ }
212+
213+
135214/*----------------------------------------------------------------------------*/
136215/* Open and allocate a mesh in serial, then read and write it in parallel */
137216/*----------------------------------------------------------------------------*/
138217
139218int main (int ArgCnt , char * * ArgVec )
140219{
141- int i , BegIdx , EndIdx ;
220+ int i , BegIdx , EndIdx , TypTab [ 3 ] ;
142221 double timer ;
143222 MshSct msh ;
144223
145224
146225 // Read the number of threads and the filenames from the command line
147- if (ArgCnt == 4 )
226+ if (ArgCnt == 5 )
148227 {
149228 msh .InpNam = (char * )* ++ ArgVec ;
150229 msh .OutNam = (char * )* ++ ArgVec ;
230+ msh .SolNam = (char * )* ++ ArgVec ;
151231 msh .NmbCpu = atoi (* ++ ArgVec );
152232 }
153233 else
154234 {
155- puts ("parallel_io InputMesh OutpuMesh NmbThreads" );
235+ puts ("parallel_io InputMesh OutputMesh OutputSolution NmbThreads" );
156236 exit (0 );
157237 }
158238
@@ -169,41 +249,46 @@ int main(int ArgCnt, char **ArgVec)
169249
170250 printf ("InpMsh : idx = %lld, version = %d, dimension = %d\n" , msh .InpMsh , msh .ver , msh .dim );
171251
172- if (msh .dim != 3 || msh . ver < 2 )
252+ if (msh .dim != 3 )
173253 {
174- puts ("Dimension must 3 and file version must be >= 2 " );
254+ puts ("Dimension must 3" );
175255 exit (1 );
176256 }
177257
178- // Read the number of vertices and tets
258+ // Read the number of vertices, triangles and tets
179259 msh .NmbVer = GmfStatKwd (msh .InpMsh , GmfVertices );
180- printf ("InpMsh : nmb vertices = %lld\n" , msh .NmbVer );
260+ printf ("InpMsh : nmb vertices = %lld\n" , msh .NmbVer );
261+
262+ msh .NmbTri = GmfStatKwd (msh .InpMsh , GmfTriangles );
263+ printf ("InpMsh : nmb triangles = %lld\n\n" , msh .NmbTri );
181264
182265 msh .NmbTet = GmfStatKwd (msh .InpMsh , GmfTetrahedra );
183- printf ("InpMsh : nmb tets = %lld\n\n" , msh .NmbTet );
266+ printf ("InpMsh : nmb tets = %lld\n\n" , msh .NmbTet );
184267
185- if (!msh .NmbVer || !msh .NmbTet )
268+ if (!msh .NmbVer || !msh .NmbTri || ! msh . NmbTet )
186269 {
187- puts ("This example only works on meshes made of vertices and tetrahedra" );
270+ puts ("This example only works on meshes made of vertices, triangles and tetrahedra" );
188271 exit (1 );
189272 }
190273
191274 // Allocate all tables
192275 msh .VerTab = malloc ((msh .NmbVer + 1 ) * 3 * sizeof (double ));
193276 msh .RefTab = malloc ((msh .NmbVer + 1 ) * sizeof (int ));
277+ msh .TriTab = malloc ((msh .NmbTri + 1 ) * 4 * sizeof (int ));
194278 msh .TetTab = malloc ((msh .NmbTet + 1 ) * 5 * sizeof (int ));
195279
196- if (!msh .VerTab || !msh .RefTab || !msh .TetTab )
280+ if (!msh .VerTab || !msh .RefTab || !msh .TriTab || ! msh . TetTab )
197281 {
198282 puts ("Failed to allocate memory" );
199283 exit (1 );
200284 }
201285
202286 msh .ParIdx = InitParallel (msh .NmbCpu );
203287 msh .VerTyp = NewType (msh .ParIdx , msh .NmbVer );
288+ msh .TriTyp = NewType (msh .ParIdx , msh .NmbTri );
204289 msh .TetTyp = NewType (msh .ParIdx , msh .NmbTet );
205290
206- if (!msh .ParIdx || !msh .VerTyp || !msh .TetTyp )
291+ if (!msh .ParIdx || !msh .VerTyp || !msh .TriTyp || ! msh . TetTyp )
207292 {
208293 puts ("Failed to initialize the parallelism with the LPlib" );
209294 exit (1 );
@@ -214,6 +299,7 @@ int main(int ArgCnt, char **ArgVec)
214299
215300 timer = GetWallClock ();
216301 LaunchParallel (msh .ParIdx , msh .VerTyp , 0 , (void * )ScaVer , (void * )& msh );
302+ LaunchParallel (msh .ParIdx , msh .TriTyp , 0 , (void * )ScaTri , (void * )& msh );
217303 LaunchParallel (msh .ParIdx , msh .TetTyp , 0 , (void * )ScaTet , (void * )& msh );
218304 printf ("Time for reading: %g seconds\n\n" , GetWallClock () - timer );
219305
@@ -242,7 +328,18 @@ int main(int ArgCnt, char **ArgVec)
242328
243329 LaunchParallel (msh .ParIdx , msh .VerTyp , 0 , (void * )RecVer , (void * )& msh );
244330
245- // Write the Tetrahedra
331+ // Write the triangles
332+ if (!(msh .OutMsh = GmfOpenMesh (msh .OutNam , GmfStopParallelWrite , msh .ver , msh .dim )))
333+ {
334+ printf ("Unable to reopen the output mesh: %s\n" , msh .OutNam );
335+ exit (1 );
336+ }
337+
338+ GmfSetKwd (msh .OutMsh , GmfTriangles , msh .NmbTri );
339+ GmfCloseUnfinishedMesh (msh .OutMsh );
340+ LaunchParallel (msh .ParIdx , msh .TriTyp , 0 , (void * )RecTri , (void * )& msh );
341+
342+ // Write the tetrahedra
246343 if (!(msh .OutMsh = GmfOpenMesh (msh .OutNam , GmfStopParallelWrite , msh .ver , msh .dim )))
247344 {
248345 printf ("Unable to reopen the output mesh: %s\n" , msh .OutNam );
@@ -253,6 +350,7 @@ int main(int ArgCnt, char **ArgVec)
253350 GmfCloseUnfinishedMesh (msh .OutMsh );
254351 LaunchParallel (msh .ParIdx , msh .TetTyp , 0 , (void * )RecTet , (void * )& msh );
255352
353+ // Final reopen in sequential mode and proposer file closing
256354 if (!(msh .OutMsh = GmfOpenMesh (msh .OutNam , GmfStopParallelWrite , msh .ver , msh .dim )))
257355 {
258356 printf ("Unable to reopen the output mesh: %s\n" , msh .OutNam );
@@ -261,6 +359,40 @@ int main(int ArgCnt, char **ArgVec)
261359
262360 GmfCloseMesh (msh .OutMsh );
263361
362+
363+ /*-----------------------------------*/
364+ /* Write the solb file */
365+ /*-----------------------------------*/
366+
367+ // Create the mesh file
368+ if (!(msh .OutSol = GmfOpenMesh (msh .SolNam , GmfWrite , msh .ver , msh .dim )))
369+ {
370+ printf ("Unable to create the output mesh: %s\n" , msh .SolNam );
371+ exit (1 );
372+ }
373+
374+ // Write the solutions at vertices
375+ TypTab [0 ] = GmfVec ;
376+ GmfSetKwd (msh .OutSol , GmfSolAtVertices , msh .NmbVer , 1 , TypTab );
377+
378+ if (!GmfCloseUnfinishedMesh (msh .OutSol ))
379+ {
380+ puts ("Parallel I/O are not available: please recompile with -DWITH_GMF_AIO" );
381+ exit (1 );
382+ }
383+
384+ LaunchParallel (msh .ParIdx , msh .VerTyp , 0 , (void * )RecSol , (void * )& msh );
385+
386+ // Final reopen in sequential mode and proposer file closing
387+ if (!(msh .OutSol = GmfOpenMesh (msh .OutNam , GmfStopParallelWrite , msh .ver , msh .dim )))
388+ {
389+ printf ("Unable to reopen the output mesh: %s\n" , msh .OutNam );
390+ exit (1 );
391+ }
392+
393+ GmfCloseMesh (msh .OutSol );
394+
395+
264396 printf ("Time for writing: %g seconds\n\n" , GetWallClock () - timer );
265397
266398 StopParallel (msh .ParIdx );
0 commit comments