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 mmgslib (basic use)
5
6
7PROGRAM main
8 IMPLICIT NONE
9
10 !> Include the mmgs library hader file
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
21 !> To save final mesh in a file
22 INTEGER :: inm=10
23 !> To manually recover the mesh
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)"
30 !> to cast integers into MMG5F_INT integers
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
46 !> ------------------------------ STEP I --------------------------
47 !! 1) Initialisation of mesh and sol structures
48 !! args of InitMesh:
49 !! MMG5_ARG_start: we start to give the args of a variadic func
50 !! MMG5_ARG_ppMesh: next arg will be a pointer over a MMG5_pMesh
51 !! mmgMesh: your MMG5_pMesh (that store your mesh)
52 !! MMG5_ARG_ppMet: next arg will be a pointer over a MMG5_pSol storing a metric
53 !! mmgSol: your MMG5_pSol (that store your metric) */
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
62 !> 2) Build mesh in MMG5 format
63 !! Two solutions: just use the MMGS_loadMesh function that will read a .mesh(b)
64 !! file formatted or manually set your mesh using the MMGS_Set* functions
65
66 !> Manually set of the mesh
67 !! a) give the size of the mesh: 12 vertices, 20 triangles, 0 edges
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
74 !> b) give the vertices: for each vertex, give the coordinates, the reference
75 !! and the position in mesh of the vertex
76 !! Note that coordinates must be in double precision to match with the coordinate
77 !! size in the C-library
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
104 !> d) give the triangles (not mandatory): for each triangle,
105 !! give the vertices index, the reference and the position of the triangle
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
147 !> 3) Build sol in MMG5 format
148 !! Two solutions: just use the MMGS_loadSol function that will read a .sol(b)
149 !! file formatted or manually set your sol using the MMGS_Set* functions
150
151 !> Manually set of the sol
152 !! a) give info for the sol structure: sol applied on vertex entities,
153 !! number of vertices=12, the sol is scalar
154 CALL mmgs_set_solsize(mmgmesh,mmgsol,mmg5_vertex,np,mmg5_scalar,ier)
155 IF ( ier /= 1 ) CALL exit(105)
156
157 !> b) give solutions values and positions
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
163 !> 4) (not mandatory): check if the number of given entities match with mesh size
164 CALL mmgs_chk_meshdata(mmgmesh,mmgsol,ier)
165 IF ( ier /= 1 ) CALL exit(107)
166
167 !> ------------------------------ STEP II --------------------------
168 !! remesh function
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
178 !> ------------------------------ STEP III --------------------------
179 !! get results */
180 !! Two solutions: just use the MMGS_saveMesh/MMGS_saveSol functions
181 !! that will write .mesh(b)/.sol formatted files or manually get your mesh/sol
182 !! using the MMGS_getMesh/MMGS_getSol functions
183
184 !> 1) Manually get the mesh
185 OPEN(unit=inm,file=trim(adjustl(fileout))//".mesh",form="formatted",status="replace")
186 WRITE(inm,*) "MeshVersionFormatted 2"
187 WRITE(inm,*) "Dimension 3"
188
189 !> a) get the new size of the mesh: vertices, triangles, edges (overwriting of old var)
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
208 !> b) Vertex recovering
209 !! Note that coordinates must be in double precision to match with the coordinate
210 !! size in the C-library
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
243 !> d) Triangles recovering
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
262 !> e) Edges recovering
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
291 !> 2) Manually get the solution
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
298 !> a) get the size of the sol: type of entity (SolAtVertices,...),
299 !! number of sol, type of solution (scalar, tensor...)
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
312 !> b) Vertex recovering
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
322 !> 3) Free the MMGS5 structures
323 CALL mmgs_free_all(mmg5_arg_start, &
324 mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
325 mmg5_arg_end)
326END 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