Mmg
Simplicial remeshers (mesh adaptation, isovalue discretization, lagrangian movement)
main.F90
Go to the documentation of this file.
1!!> @author
2!> Cecile Dobrzynski, Charles Dapogny, Pascal Frey and Algiane Froehly
3!> @brief
4!> Example for using mmg3dlib (basic use)
5
6PROGRAM main
7
8 IMPLICIT NONE
9
10!> Include the mmg3d library hader file
11! if the header file is in the "include" directory
12! #include "libmmg3df.h"
13! if the header file is in "include/mmg/mmg3d"
14#include "mmg/mmg3d/libmmg3df.h"
15
16 mmg5_data_ptr_t :: mmgmesh
17 mmg5_data_ptr_t :: mmgsol
18 INTEGER :: ier,argc
19 CHARACTER(len=300) :: exec_name,fileout
20
21 !> To save final mesh in a file
22 INTEGER :: inm=10
23 !> To manually recover the mesh
24 INTEGER(MMG5F_INT):: k, np, ne, nt, na, nc, nr, nreq, ref
25 INTEGER(MMG5F_INT):: tetra(4), tria(3), edge(2)
26 INTEGER :: typentity, typsol
27 DOUBLE PRECISION :: point(3),sol
28 INTEGER, DIMENSION(:), ALLOCATABLE :: corner, required, ridge
29 CHARACTER(LEN=31) :: fmt="(E14.8,1X,E14.8,1X,E14.8,1X,I3)"
30 INTEGER(MMG5F_INT),DIMENSION(2) :: ktet
31 INTEGER,DIMENSION(2) :: iface
32 !> to cast integers into MMG5F_INT integers
33 INTEGER,PARAMETER :: immg = mmg5f_int
34
35 print*," -- TEST MMG3DLIB"
36
37 argc = command_argument_count();
38 CALL get_command_argument(0, exec_name)
39
40 IF ( argc /=1 ) THEN
41 print*," Usage: ",trim(exec_name)," output_filename"
42 CALL exit(1);
43 ENDIF
44
45 ! Name and path of the mesh file
46 CALL get_command_argument(1, fileout)
47
48 !> ------------------------------ STEP I --------------------------
49 !! 1) Initialisation of mesh and sol structures
50 !! args of InitMesh:
51 !! MMG5_ARG_start: we start to give the args of a variadic func
52 !! MMG5_ARG_ppMesh: next arg will be a pointer over a MMG5_pMesh
53 !! mmgMesh: your MMG5_pMesh (that store your mesh)
54 !! MMG5_ARG_ppMet: next arg will be a pointer over a MMG5_pSol storing a metric
55 !! mmgSol: your MMG5_pSol (that store your metric) */
56
57 mmgmesh = 0
58 mmgsol = 0
59
60 CALL mmg3d_init_mesh(mmg5_arg_start, &
61 mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
62 mmg5_arg_end)
63 CALL mmg3d_set_iparameter(mmgmesh,mmgsol,mmg3d_iparam_verbose,5_immg,ier)
64
65 !> 2) Build mesh in MMG5 format
66 !! Two solutions: just use the MMG3D_loadMesh function that will read a .mesh(b)
67 !! file formatted or manually set your mesh using the MMG3D_Set* functions
68
69 !> Manually set of the mesh
70 !! a) give the size of the mesh: 12 vertices, 12 tetra,0 prisms, 20 triangles,
71 !! 0 quads, 0 edges
72 np = 12
73 ne = 12
74 nt = 20
75 CALL mmg3d_set_meshsize(mmgmesh,np,ne,0_immg,nt,0_immg,0_immg,ier)
76 IF ( ier /= 1 ) CALL exit(101)
77
78 !> b) give the vertices: for each vertex, give the coordinates, the reference
79 !! and the position in mesh of the vertex
80 !! Note that coordinates must be in double precision to match with the coordinate
81 !! size in the C-library
82
83 CALL mmg3d_set_vertex(mmgmesh, 0.0d0, 0.0d0, 0.0d0, 0_immg, 1_immg,ier)
84 IF ( ier /= 1 ) CALL exit(102)
85 CALL mmg3d_set_vertex(mmgmesh, 0.5d0, 0.0d0, 0.0d0, 0_immg, 2_immg,ier)
86 IF ( ier /= 1 ) CALL exit(102)
87 CALL mmg3d_set_vertex(mmgmesh, 0.5d0, 0.0d0, 1.0d0, 0_immg, 3_immg,ier)
88 IF ( ier /= 1 ) CALL exit(102)
89 CALL mmg3d_set_vertex(mmgmesh, 0.0d0, 0.0d0, 1.0d0, 0_immg, 4_immg,ier)
90 IF ( ier /= 1 ) CALL exit(102)
91 CALL mmg3d_set_vertex(mmgmesh, 0.0d0, 1.0d0, 0.0d0, 0_immg, 5_immg,ier)
92 IF ( ier /= 1 ) CALL exit(102)
93 CALL mmg3d_set_vertex(mmgmesh, 0.5d0, 1.0d0, 0.0d0, 0_immg, 6_immg,ier)
94 IF ( ier /= 1 ) CALL exit(102)
95 CALL mmg3d_set_vertex(mmgmesh, 0.5d0, 1.0d0, 1.0d0, 0_immg, 7_immg,ier)
96 IF ( ier /= 1 ) CALL exit(102)
97 CALL mmg3d_set_vertex(mmgmesh, 0.0d0, 1.0d0, 1.0d0, 0_immg, 8_immg,ier)
98 IF ( ier /= 1 ) CALL exit(102)
99 CALL mmg3d_set_vertex(mmgmesh, 1.0d0, 0.0d0, 0.0d0, 0_immg, 9_immg,ier)
100 IF ( ier /= 1 ) CALL exit(102)
101 CALL mmg3d_set_vertex(mmgmesh, 1.0d0, 1.0d0, 0.0d0, 0_immg, 10_immg,ier)
102 IF ( ier /= 1 ) CALL exit(102)
103 CALL mmg3d_set_vertex(mmgmesh, 1.0d0, 0.0d0, 1.0d0, 0_immg, 11_immg,ier)
104 IF ( ier /= 1 ) CALL exit(102)
105 CALL mmg3d_set_vertex(mmgmesh, 1.0d0, 1.0d0, 1.0d0, 0_immg, 12_immg,ier)
106 IF ( ier /= 1 ) CALL exit(102)
107
108 !> c) give the tetrahedras: for each tetrahedra,
109 !! give the vertices index, the reference and the position of the tetra
110 ref = 1
111 CALL mmg3d_set_tetrahedron(mmgmesh, 1_immg, 4_immg, 2_immg, 8_immg,ref, 1_immg,ier)
112 IF ( ier /= 1 ) CALL exit(103)
113 CALL mmg3d_set_tetrahedron(mmgmesh, 8_immg, 3_immg, 2_immg, 7_immg,ref, 2_immg,ier)
114 IF ( ier /= 1 ) CALL exit(103)
115 CALL mmg3d_set_tetrahedron(mmgmesh, 5_immg, 2_immg, 6_immg, 8_immg,ref, 3_immg,ier)
116 IF ( ier /= 1 ) CALL exit(103)
117 CALL mmg3d_set_tetrahedron(mmgmesh, 5_immg, 8_immg, 1_immg, 2_immg,ref, 4_immg,ier)
118 IF ( ier /= 1 ) CALL exit(103)
119 CALL mmg3d_set_tetrahedron(mmgmesh, 7_immg, 2_immg, 8_immg, 6_immg,ref, 5_immg,ier)
120 IF ( ier /= 1 ) CALL exit(103)
121 CALL mmg3d_set_tetrahedron(mmgmesh, 2_immg, 4_immg, 3_immg, 8_immg,ref, 6_immg,ier)
122 IF ( ier /= 1 ) CALL exit(103)
123 ref = 2
124 CALL mmg3d_set_tetrahedron(mmgmesh, 9_immg, 2_immg, 3_immg, 7_immg,ref, 7_immg,ier)
125 IF ( ier /= 1 ) CALL exit(103)
126 CALL mmg3d_set_tetrahedron(mmgmesh, 7_immg, 11_immg, 9_immg, 12_immg,ref, 8_immg,ier)
127 IF ( ier /= 1 ) CALL exit(103)
128 CALL mmg3d_set_tetrahedron(mmgmesh, 6_immg, 9_immg,10_immg, 7_immg,ref, 9_immg,ier)
129 IF ( ier /= 1 ) CALL exit(103)
130 CALL mmg3d_set_tetrahedron(mmgmesh, 6_immg, 7_immg, 2_immg, 9_immg,ref,10_immg,ier)
131 IF ( ier /= 1 ) CALL exit(103)
132 CALL mmg3d_set_tetrahedron(mmgmesh, 12_immg, 9_immg, 7_immg, 10_immg,ref,11_immg,ier)
133 IF ( ier /= 1 ) CALL exit(103)
134 CALL mmg3d_set_tetrahedron(mmgmesh, 9_immg, 3_immg,11_immg, 7_immg,ref,12_immg,ier)
135 IF ( ier /= 1 ) CALL exit(103)
136
137 !> d) give the triangles (not mandatory): for each triangle,
138 !! give the vertices index, the reference and the position of the triangle
139 ref = 3
140 CALL mmg3d_set_triangle(mmgmesh, 1_immg, 4_immg, 8_immg, ref, 1_immg,ier)
141 IF ( ier /= 1 ) CALL exit(104)
142 CALL mmg3d_set_triangle(mmgmesh, 1_immg, 2_immg, 4_immg, ref, 2_immg,ier)
143 IF ( ier /= 1 ) CALL exit(104)
144 CALL mmg3d_set_triangle(mmgmesh, 8_immg, 3_immg, 7_immg, ref, 3_immg,ier)
145 IF ( ier /= 1 ) CALL exit(104)
146 CALL mmg3d_set_triangle(mmgmesh, 5_immg, 8_immg, 6_immg, ref, 4_immg,ier)
147 IF ( ier /= 1 ) CALL exit(104)
148 CALL mmg3d_set_triangle(mmgmesh, 5_immg, 6_immg, 2_immg, ref, 5_immg,ier)
149 IF ( ier /= 1 ) CALL exit(104)
150 CALL mmg3d_set_triangle(mmgmesh, 5_immg, 2_immg, 1_immg, ref, 6_immg,ier)
151 IF ( ier /= 1 ) CALL exit(104)
152 CALL mmg3d_set_triangle(mmgmesh, 5_immg, 1_immg, 8_immg, ref, 7_immg,ier)
153 IF ( ier /= 1 ) CALL exit(104)
154 CALL mmg3d_set_triangle(mmgmesh, 7_immg, 6_immg, 8_immg, ref, 8_immg,ier)
155 IF ( ier /= 1 ) CALL exit(104)
156 CALL mmg3d_set_triangle(mmgmesh, 4_immg, 3_immg, 8_immg, ref, 9_immg,ier)
157 IF ( ier /= 1 ) CALL exit(104)
158 CALL mmg3d_set_triangle(mmgmesh, 2_immg, 3_immg, 4_immg, ref,10_immg,ier)
159 IF ( ier /= 1 ) CALL exit(104)
160 ref = 4
161 CALL mmg3d_set_triangle(mmgmesh, 9_immg, 3_immg, 2_immg, ref,11_immg,ier)
162 IF ( ier /= 1 ) CALL exit(104)
163 CALL mmg3d_set_triangle(mmgmesh, 11_immg, 9_immg,12_immg, ref,12_immg,ier)
164 IF ( ier /= 1 ) CALL exit(104)
165 CALL mmg3d_set_triangle(mmgmesh, 7_immg, 11_immg,12_immg, ref,13_immg,ier)
166 IF ( ier /= 1 ) CALL exit(104)
167 CALL mmg3d_set_triangle(mmgmesh, 6_immg, 7_immg,10_immg, ref,14_immg,ier)
168 IF ( ier /= 1 ) CALL exit(104)
169 CALL mmg3d_set_triangle(mmgmesh, 6_immg, 10_immg, 9_immg, ref,15_immg,ier)
170 IF ( ier /= 1 ) CALL exit(104)
171 CALL mmg3d_set_triangle(mmgmesh, 6_immg, 9_immg, 2_immg, ref,16_immg,ier)
172 IF ( ier /= 1 ) CALL exit(104)
173 CALL mmg3d_set_triangle(mmgmesh, 12_immg, 10_immg, 7_immg, ref,17_immg,ier)
174 IF ( ier /= 1 ) CALL exit(104)
175 CALL mmg3d_set_triangle(mmgmesh, 12_immg, 9_immg,10_immg, ref,18_immg,ier)
176 IF ( ier /= 1 ) CALL exit(104)
177 CALL mmg3d_set_triangle(mmgmesh, 3_immg, 11_immg, 7_immg, ref,19_immg,ier)
178 IF ( ier /= 1 ) CALL exit(104)
179 CALL mmg3d_set_triangle(mmgmesh, 9_immg, 11_immg, 3_immg, ref,20_immg,ier)
180 IF ( ier /= 1 ) CALL exit(104)
181
182 !> 3) Build sol in MMG5 format
183 !! Two solutions: just use the MMG3D_loadSol function that will read a .sol(b)
184 !! file formatted or manually set your sol using the MMG3D_Set* functions
185
186 !> Manually set of the sol
187 !! a) give info for the sol structure: sol applied on vertex entities,
188 !! number of vertices=12, the sol is scalar
189 np = 12
190 CALL mmg3d_set_solsize(mmgmesh,mmgsol,mmg5_vertex,np,mmg5_scalar,ier)
191 IF ( ier /= 1 ) CALL exit(105)
192
193 !> b) give solutions values and positions
194 DO k=1,12
195 CALL mmg3d_set_scalarsol(mmgsol,0.5d0,k,ier)
196 IF ( ier /= 1 ) CALL exit(106)
197 ENDDO
198
199 !> 4) (not mandatory): check if the number of given entities match with mesh size
200 CALL mmg3d_chk_meshdata(mmgmesh,mmgsol,ier)
201 IF ( ier /= 1 ) CALL exit(107)
202
203 !> ------------------------------ STEP II --------------------------
204 !! remesh function
205 CALL mmg3d_mmg3dlib(mmgmesh,mmgsol,ier)
206
207 IF ( ier == mmg5_strongfailure ) THEN
208 print*,"BAD ENDING OF MMG3DLIB: UNABLE TO SAVE MESH"
209 stop mmg5_strongfailure
210 ELSE IF ( ier == mmg5_lowfailure ) THEN
211 print*,"BAD ENDING OF MMG3DLIB"
212 ENDIF
213
214 !> ------------------------------ STEP III --------------------------
215 !! get results */
216 !! Two solutions: just use the MMG3D_saveMesh/MMG3D_saveSol functions
217 !! that will write .mesh(b)/.sol formatted files or manually get your mesh/sol
218 !! using the MMG3D_getMesh/MMG3D_getSol functions
219
220 !> 1) Manually get the mesh
221 OPEN(unit=inm,file=trim(adjustl(fileout))//".mesh",form="formatted",status="replace")
222 WRITE(inm,*) "MeshVersionFormatted 2"
223 WRITE(inm,*) "Dimension 3"
224
225 !> a) get the size of the mesh: vertices, tetra,prisms, triangles, quads,edges
226 CALL mmg3d_get_meshsize(mmgmesh,np,ne,%val(0_immg),nt,%val(0_immg),na,ier)
227 IF ( ier /= 1 ) CALL exit(108)
228
229 ! Table to know if a vertex is corner
230 ALLOCATE(corner(np))
231
232 ! Table to know if a vertex/tetra/tria/edge is required
233 ALLOCATE(required(max(max(np,ne),max(nt,na))))
234
235 ! Table to know if a coponant is corner and/or required
236 ALLOCATE(ridge(na))
237
238 nreq = 0; nc = 0
239 WRITE(inm,*)
240 WRITE(inm,*) "Vertices"
241 WRITE(inm,*) np
242
243 DO k=1, np
244 !> b) Vertex recovering
245 !! Note that coordinates must be in double precision to match with the coordinate
246 !! size in the C-library
247 CALL mmg3d_get_vertex(mmgmesh,point(1),point(2),point(3),&
248 ref,corner(k),required(k),ier)
249 IF ( ier /= 1 ) CALL exit(109)
250
251 WRITE(inm,fmt) point(1),point(2),point(3),ref
252 IF ( corner(k)/=0 ) nc=nc+1
253 IF ( required(k)/=0 ) nreq=nreq+1
254 ENDDO
255
256 WRITE(inm,*)
257 WRITE(inm,*) "Corners"
258 WRITE(inm,*) nc
259
260 DO k=1, np
261 IF ( corner(k)/=0 ) WRITE(inm,*) k
262 ENDDO
263 WRITE(inm,*)
264
265 WRITE(inm,*) "RequiredVertices"
266 WRITE(inm,*) nreq
267
268 DO k=1,np
269 IF ( required(k)/=0 ) WRITE(inm,*) k
270 ENDDO
271 WRITE(inm,*)
272 DEALLOCATE(corner)
273
274 nreq = 0;
275 WRITE(inm,*) "Triangles"
276 WRITE(inm,*) nt
277
278 DO k=1,nt
279 !> d) Triangles recovering
280 CALL mmg3d_get_triangle(mmgmesh,tria(1),tria(2),tria(3),ref,required(k),ier)
281 IF ( ier /= 1 ) CALL exit(110)
282 WRITE(inm,*) tria(1),tria(2),tria(3),ref
283 IF ( required(k)/=0 ) nreq=nreq+1;
284 ENDDO
285 WRITE(inm,*)
286
287 WRITE(inm,*) "RequiredTriangles"
288 WRITE(inm,*) nreq
289 DO k=1,nt
290 IF ( required(k)/=0 ) WRITE(inm,*) k
291 ENDDO
292 WRITE(inm,*)
293
294 nreq = 0;nr = 0;
295 WRITE(inm,*) "Edges"
296 WRITE(inm,*) na
297 DO k=1,na
298 !> e) Edges recovering
299 CALL mmg3d_get_edge(mmgmesh,edge(1),edge(2),ref,ridge(k),required(k),ier)
300 IF ( ier /= 1 ) CALL exit(111)
301 WRITE(inm,*) edge(1),edge(2),ref
302 IF ( ridge(k)/=0 ) nr = nr+1
303 IF ( required(k)/=0 ) nreq = nreq+1
304 ENDDO
305 WRITE(inm,*)
306
307 WRITE(inm,*) "RequiredEdges"
308 WRITE(inm,*) nreq
309 DO k=1,na
310 IF ( required(k) /=0 ) WRITE(inm,*) k
311 ENDDO
312 WRITE(inm,*)
313
314 WRITE(inm,*) "Ridges"
315 WRITE(inm,*) nr
316 DO k=1,na
317 IF ( ridge(k) /=0 ) WRITE(inm,*) k
318 ENDDO
319 WRITE(inm,*)
320
321 nreq = 0;
322 WRITE(inm,*) "Tetrahedra"
323 WRITE(inm,*) ne
324 DO k=1,ne
325 !> c) Tetra recovering
326 CALL mmg3d_get_tetrahedron(mmgmesh,tetra(1),tetra(2),tetra(3),tetra(4),&
327 ref,required(k),ier)
328 IF ( ier /= 1 ) CALL exit(112)
329 WRITE(inm,*) tetra(1),tetra(2),tetra(3),tetra(4),ref
330 IF ( required(k) /= 0 ) nreq = nreq+1
331 ENDDO
332 WRITE(inm,*)
333
334 WRITE(inm,*) "RequiredTetrahedra"
335 WRITE(inm,*) nreq
336 DO k=1,ne
337 IF ( required(k) /= 0 ) WRITE(inm,*) k
338 ENDDO
339
340 WRITE(inm,*) "End"
341 CLOSE(inm)
342
343 DEALLOCATE(required)
344 DEALLOCATE(ridge)
345
346 !> 2) Manually get the solution
347 OPEN(unit=inm,file=trim(adjustl(fileout))//".sol",form="formatted",status="replace")
348 WRITE(inm,*) "MeshVersionFormatted 2"
349 WRITE(inm,*) "Dimension 3"
350 WRITE(inm,*)
351
352 !> a) get the size of the sol: type of entity (SolAtVertices,...),
353 !! number of sol, type of solution (scalar, tensor...)
354 CALL mmg3d_get_solsize(mmgmesh,mmgsol,typentity,np,typsol,ier)
355 IF ( ier /= 1 ) CALL exit(113)
356
357 IF ( ( typentity /= mmg5_vertex ) .OR. ( typsol /= mmg5_scalar ) ) THEN
358 CALL exit(114);
359 ENDIF
360
361 WRITE(inm,*) "SolAtVertices"
362 WRITE(inm,*) np
363 WRITE(inm,*) "1 1"
364 WRITE(inm,*)
365 DO k=1,np
366 !> b) Vertex recovering
367 CALL mmg3d_get_scalarsol(mmgsol,sol,ier)
368 IF ( ier /= 1 ) CALL exit(115)
369 WRITE(inm,*) sol
370 ENDDO
371 WRITE(inm,*)
372
373 WRITE(inm,*) "End"
374 CLOSE(inm)
375
376 !> 3) Free the MMG3D5 structures
377 CALL mmg3d_free_all(mmg5_arg_start, &
378 mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
379 mmg5_arg_end)
380END PROGRAM main
main
int main(int argc, char *argv[])
Definition mmg2d.c:275
Generated on Mon Apr 7 2025 12:27:12 for Mmg by doxygen 1.13.2