Mmg
Simplicial remeshers (mesh adaptation, isovalue discretization, lagrangian movement)
main.F90
Go to the documentation of this file.
1!!> @author
5
6PROGRAM main
7
8 IMPLICIT NONE
9
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
22 INTEGER :: inm=10
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
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
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
68
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
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
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
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
185
189 np = 12
190 CALL mmg3d_set_solsize(mmgmesh,mmgsol,mmg5_vertex,np,mmg5_scalar,ier)
191 IF ( ier /= 1 ) CALL exit(105)
192
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
200 CALL mmg3d_chk_meshdata(mmgmesh,mmgsol,ier)
201 IF ( ier /= 1 ) CALL exit(107)
202
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
219
221 OPEN(unit=inm,file=trim(adjustl(fileout))//".mesh",form="formatted",status="replace")
222 WRITE(inm,*) "MeshVersionFormatted 2"
223 WRITE(inm,*) "Dimension 3"
224
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
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
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
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
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
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
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
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
377 CALL mmg3d_free_all(mmg5_arg_start, &
378 mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
379 mmg5_arg_end)
380END PROGRAM main
int main(int argc, char *argv[])
Definition: mmg2d.c:275