Mmg
Simplicial remeshers (mesh adaptation, isovalue discretization, lagrangian movement)
main.F90
Go to the documentation of this file.
1!!> @author
5
6
7PROGRAM main
8 IMPLICIT NONE
9
11 ! if the header file is in the "include" directory
12 ! #include "libmmgsf.h"
13 ! if the header file is in "include/mmg/mmgs"
14#include "mmg/mmgs/libmmgsf.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, nt, na, nc, nr, nreq
25 INTEGER :: typentity, typsol
26 INTEGER(MMG5F_INT):: ref, tria(3), edge(2)
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)"
31 INTEGER,PARAMETER :: immg = mmg5f_int
32
33 print*," -- TEST MMGSLIB"
34
35 argc = command_argument_count();
36 CALL get_command_argument(0, exec_name)
37
38 IF ( argc /=1 ) THEN
39 print*," Usage: ",trim(exec_name)," output_file_name"
40 CALL exit(1);
41 ENDIF
42
43 ! Name and path of the mesh file
44 CALL get_command_argument(1, fileout)
45
54
55 mmgmesh = 0
56 mmgsol = 0
57
58 CALL mmgs_init_mesh(mmg5_arg_start, &
59 mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
60 mmg5_arg_end)
61
65
68 np = 12
69 nt = 20
70 na = 0
71 CALL mmgs_set_meshsize(mmgmesh,np,nt,na,ier)
72 IF ( ier /= 1 ) CALL exit(101)
73
78 ref = 0
79 CALL mmgs_set_vertex(mmgmesh, 0.0d0, 0.0d0, 0.0d0, ref, 1_immg,ier)
80 IF ( ier /= 1 ) CALL exit(102)
81 CALL mmgs_set_vertex(mmgmesh, 0.5d0, 0.0d0, 0.0d0, ref, 2_immg,ier)
82 IF ( ier /= 1 ) CALL exit(102)
83 CALL mmgs_set_vertex(mmgmesh, 0.5d0, 0.0d0, 1.0d0, ref, 3_immg,ier)
84 IF ( ier /= 1 ) CALL exit(102)
85 CALL mmgs_set_vertex(mmgmesh, 0.0d0, 0.0d0, 1.0d0, ref, 4_immg,ier)
86 IF ( ier /= 1 ) CALL exit(102)
87 CALL mmgs_set_vertex(mmgmesh, 0.0d0, 1.0d0, 0.0d0, ref, 5_immg,ier)
88 IF ( ier /= 1 ) CALL exit(102)
89 CALL mmgs_set_vertex(mmgmesh, 0.5d0, 1.0d0, 0.0d0, ref, 6_immg,ier)
90 IF ( ier /= 1 ) CALL exit(102)
91 CALL mmgs_set_vertex(mmgmesh, 0.5d0, 1.0d0, 1.0d0, ref, 7_immg,ier)
92 IF ( ier /= 1 ) CALL exit(102)
93 CALL mmgs_set_vertex(mmgmesh, 0.0d0, 1.0d0, 1.0d0, ref, 8_immg,ier)
94 IF ( ier /= 1 ) CALL exit(102)
95 CALL mmgs_set_vertex(mmgmesh, 1.0d0, 0.0d0, 0.0d0, ref, 9_immg,ier)
96 IF ( ier /= 1 ) CALL exit(102)
97 CALL mmgs_set_vertex(mmgmesh, 1.0d0, 1.0d0, 0.0d0, ref, 10_immg,ier)
98 IF ( ier /= 1 ) CALL exit(102)
99 CALL mmgs_set_vertex(mmgmesh, 1.0d0, 0.0d0, 1.0d0, ref, 11_immg,ier)
100 IF ( ier /= 1 ) CALL exit(102)
101 CALL mmgs_set_vertex(mmgmesh, 1.0d0, 1.0d0, 1.0d0, ref, 12_immg,ier)
102 IF ( ier /= 1 ) CALL exit(102)
103
106 CALL mmgs_set_triangle(mmgmesh, 1_immg, 4_immg, 8_immg, 3_immg, 1_immg,ier)
107 IF ( ier /= 1 ) CALL exit(104)
108 CALL mmgs_set_triangle(mmgmesh, 1_immg, 2_immg, 4_immg, 3_immg, 2_immg,ier)
109 IF ( ier /= 1 ) CALL exit(104)
110 CALL mmgs_set_triangle(mmgmesh, 8_immg, 3_immg, 7_immg, 3_immg, 3_immg,ier)
111 IF ( ier /= 1 ) CALL exit(104)
112 CALL mmgs_set_triangle(mmgmesh, 5_immg, 8_immg, 6_immg, 3_immg, 4_immg,ier)
113 IF ( ier /= 1 ) CALL exit(104)
114 CALL mmgs_set_triangle(mmgmesh, 5_immg, 6_immg, 2_immg, 3_immg, 5_immg,ier)
115 IF ( ier /= 1 ) CALL exit(104)
116 CALL mmgs_set_triangle(mmgmesh, 5_immg, 2_immg, 1_immg, 3_immg, 6_immg,ier)
117 IF ( ier /= 1 ) CALL exit(104)
118 CALL mmgs_set_triangle(mmgmesh, 5_immg, 1_immg, 8_immg, 3_immg, 7_immg,ier)
119 IF ( ier /= 1 ) CALL exit(104)
120 CALL mmgs_set_triangle(mmgmesh, 7_immg, 6_immg, 8_immg, 3_immg, 8_immg,ier)
121 IF ( ier /= 1 ) CALL exit(104)
122 CALL mmgs_set_triangle(mmgmesh, 4_immg, 3_immg, 8_immg, 3_immg, 9_immg,ier)
123 IF ( ier /= 1 ) CALL exit(104)
124 CALL mmgs_set_triangle(mmgmesh, 2_immg, 3_immg, 4_immg, 3_immg,10_immg,ier)
125 IF ( ier /= 1 ) CALL exit(104)
126 CALL mmgs_set_triangle(mmgmesh, 9_immg, 3_immg, 2_immg, 4_immg,11_immg,ier)
127 IF ( ier /= 1 ) CALL exit(104)
128 CALL mmgs_set_triangle(mmgmesh, 11_immg, 9_immg, 12_immg, 4_immg,12_immg,ier)
129 IF ( ier /= 1 ) CALL exit(104)
130 CALL mmgs_set_triangle(mmgmesh, 7_immg, 11_immg, 12_immg, 4_immg,13_immg,ier)
131 IF ( ier /= 1 ) CALL exit(104)
132 CALL mmgs_set_triangle(mmgmesh, 6_immg, 7_immg, 10_immg, 4_immg,14_immg,ier)
133 IF ( ier /= 1 ) CALL exit(104)
134 CALL mmgs_set_triangle(mmgmesh, 6_immg, 10_immg, 9_immg, 4_immg,15_immg,ier)
135 IF ( ier /= 1 ) CALL exit(104)
136 CALL mmgs_set_triangle(mmgmesh, 6_immg, 9_immg, 2_immg, 4_immg,16_immg,ier)
137 IF ( ier /= 1 ) CALL exit(104)
138 CALL mmgs_set_triangle(mmgmesh, 12_immg, 10_immg, 7_immg, 4_immg,17_immg,ier)
139 IF ( ier /= 1 ) CALL exit(104)
140 CALL mmgs_set_triangle(mmgmesh, 12_immg, 9_immg, 10_immg, 4_immg,18_immg,ier)
141 IF ( ier /= 1 ) CALL exit(104)
142 CALL mmgs_set_triangle(mmgmesh, 3_immg, 11_immg, 7_immg, 4_immg,19_immg,ier)
143 IF ( ier /= 1 ) CALL exit(104)
144 CALL mmgs_set_triangle(mmgmesh, 9_immg, 11_immg, 3_immg, 4_immg,20_immg,ier)
145 IF ( ier /= 1 ) CALL exit(104)
146
150
154 CALL mmgs_set_solsize(mmgmesh,mmgsol,mmg5_vertex,np,mmg5_scalar,ier)
155 IF ( ier /= 1 ) CALL exit(105)
156
158 DO k=1,12
159 CALL mmgs_set_scalarsol(mmgsol,0.5d0,k,ier)
160 IF ( ier /= 1 ) CALL exit(106)
161 ENDDO
162
164 CALL mmgs_chk_meshdata(mmgmesh,mmgsol,ier)
165 IF ( ier /= 1 ) CALL exit(107)
166
169 CALL mmgs_mmgslib(mmgmesh,mmgsol,ier)
170
171 IF ( ier == mmg5_strongfailure ) THEN
172 print*,"BAD ENDING OF MMGSLIB: UNABLE TO SAVE MESH"
173 stop mmg5_strongfailure
174 ELSE IF ( ier == mmg5_lowfailure ) THEN
175 print*,"BAD ENDING OF MMGSLIB"
176 ENDIF
177
183
185 OPEN(unit=inm,file=trim(adjustl(fileout))//".mesh",form="formatted",status="replace")
186 WRITE(inm,*) "MeshVersionFormatted 2"
187 WRITE(inm,*) "Dimension 3"
188
190 CALL mmgs_get_meshsize(mmgmesh,np,nt,na,ier)
191 IF ( ier /= 1 ) CALL exit(108)
192
193 ! Table to know if a vertex is corner
194 ALLOCATE(corner(np))
195
196 ! Table to know if a vertex/tetra/tria/edge is required
197 ALLOCATE(required(max(np,max(nt,na))))
198
199 ! Table to know if a coponant is corner and/or required
200 ALLOCATE(ridge(na))
201
202 nreq = 0; nc = 0
203 WRITE(inm,*)
204 WRITE(inm,*) "Vertices"
205 WRITE(inm,*) np
206
207 DO k=1, np
211 CALL mmgs_get_vertex(mmgmesh,point(1),point(2),point(3),&
212 ref,corner(k),required(k),ier)
213 IF ( ier /= 1 ) CALL exit(109)
214
215 WRITE(inm,fmt) point(1),point(2),point(3),ref
216 IF ( corner(k)/=0 ) nc=nc+1
217 IF ( required(k)/=0 ) nreq=nreq+1
218 ENDDO
219
220 WRITE(inm,*)
221 WRITE(inm,*) "Corners"
222 WRITE(inm,*) nc
223
224 DO k=1, np
225 IF ( corner(k)/=0 ) WRITE(inm,*) k
226 ENDDO
227 WRITE(inm,*)
228
229 WRITE(inm,*) "RequiredVertices"
230 WRITE(inm,*) nreq
231
232 DO k=1,np
233 IF ( required(k)/=0 ) WRITE(inm,*) k
234 ENDDO
235 WRITE(inm,*)
236 DEALLOCATE(corner)
237
238 nreq = 0;
239 WRITE(inm,*) "Triangles"
240 WRITE(inm,*) nt
241
242 DO k=1,nt
244 CALL mmgs_get_triangle(mmgmesh,tria(1),tria(2),tria(3),ref,required(k),ier)
245 IF ( ier /= 1 ) CALL exit(110)
246 WRITE(inm,*) tria(1),tria(2),tria(3),ref
247 IF ( required(k)/=0 ) nreq=nreq+1;
248 ENDDO
249 WRITE(inm,*)
250
251 WRITE(inm,*) "RequiredTriangles"
252 WRITE(inm,*) nreq
253 DO k=1,nt
254 IF ( required(k)/=0 ) WRITE(inm,*) k
255 ENDDO
256 WRITE(inm,*)
257
258 nreq = 0;nr = 0;
259 WRITE(inm,*) "Edges"
260 WRITE(inm,*) na
261 DO k=1,na
263 CALL mmgs_get_edge(mmgmesh,edge(1),edge(2),ref,ridge(k),required(k),ier)
264 IF ( ier /= 1 ) CALL exit(111)
265 WRITE(inm,*) edge(1),edge(2),ref
266 IF ( ridge(k)/=0 ) nr = nr+1
267 IF ( required(k)/=0 ) nreq = nreq+1
268 ENDDO
269 WRITE(inm,*)
270
271 WRITE(inm,*) "RequiredEdges"
272 WRITE(inm,*) nreq
273 DO k=1,na
274 IF ( required(k) /=0 ) WRITE(inm,*) k
275 ENDDO
276 WRITE(inm,*)
277
278 WRITE(inm,*) "Ridges"
279 WRITE(inm,*) nr
280 DO k=1,na
281 IF ( ridge(k) /=0 ) WRITE(inm,*) k
282 ENDDO
283 WRITE(inm,*)
284
285 WRITE(inm,*) "End"
286 CLOSE(inm)
287
288 DEALLOCATE(required)
289 DEALLOCATE(ridge)
290
292 OPEN(unit=inm,file=trim(adjustl(fileout))//".sol", &
293 & form="formatted",status="replace")
294 WRITE(inm,*) "MeshVersionFormatted 2"
295 WRITE(inm,*) "Dimension 3"
296 WRITE(inm,*)
297
300 CALL mmgs_get_solsize(mmgmesh,mmgsol,typentity,np,typsol,ier)
301 IF ( ier /= 1 ) CALL exit(113)
302
303 IF ( ( typentity /= mmg5_vertex ) .OR. ( typsol /= mmg5_scalar ) ) THEN
304 CALL exit(114);
305 ENDIF
306
307 WRITE(inm,*) "SolAtVertices"
308 WRITE(inm,*) np
309 WRITE(inm,*) "1 1"
310 WRITE(inm,*)
311 DO k=1,np
313 CALL mmgs_get_scalarsol(mmgsol,sol,ier)
314 IF ( ier /= 1 ) CALL exit(115)
315 WRITE(inm,*) sol
316 ENDDO
317 WRITE(inm,*)
318
319 WRITE(inm,*) "End"
320 CLOSE(inm)
321
323 CALL mmgs_free_all(mmg5_arg_start, &
324 mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
325 mmg5_arg_end)
326END PROGRAM main
int main(int argc, char *argv[])
Definition mmg2d.c:275