Choral
mesh_mmg.F90
Go to the documentation of this file.
1 
16 
17 module mesh_mmg
18 
21  use basic_tools
22  use real_type
23  use graph_mod
24  use cell_mod
25  use mesh_mod
26 
27  implicit none
28  private
29 
30 #ifdef WMMG
31 #include "mmg/mmg2d/libmmg2df.h"
32 
33  public :: mmg2d_retrievemesh, mmg2d_passmesh
34 #endif
35 
36  public :: remesh_levelset
37 
38 
39 contains
40 
41  !! Remesh respectively to a levelSet
42  !!
43  !! m = initial mesh / output mesh
44  !! val = level values at the vertexes of the mesh 'm'
45  !!
46  !! Output mesh has edges along the level set 0
47  !!
48  !! Optional parameters
49  !! verb = verbisity
50  !! hausd = hausdorff distance
51  !! hMax = maximal edge length
52  !!
53  subroutine remesh_levelset(m, val, verb, hausd, hMax)
54  type(mesh), intent(inout) :: m
55  real(DP) , dimension(:) :: val
56 
57  integer , optional, intent(in) :: verb
58  real(DP), optional, intent(in) :: hausd
59  real(DP), optional, intent(in) :: hMax
60 
61 #ifdef WMMG
62  integer :: ii, ier
63  mmg5_data_ptr_t :: mmgmesh, mmgsol
64 
65  if (choral_verb>0) write(*,*) &
66  & "mesh_mmg : reMesh_levelSet"
67 
68  if (.NOT.valid(m)) call quit("mesh_mmg: reMesh_levelSet:&
69  & mesh 'm' not valid")
70 
71  if (size(val,1) /= m%nbNd) call quit("mesh_mmg: reMesh_levelSet:&
72  & wrong size for 'val'")
73 
74 
75  !! Adapt the code to the mesh dimension
76  select case(m%dim)
77  case(2)
78  if ( m%cell_Count(cell_trg)==0 ) call quit(&
79  & "mesh_mmg: reMesh_levelSet:&
80  & dimension 2 mesh withour triangles")
81 
82  !! #####################
83  !!
84  !! DATA TRANSFER TO MMG
85  !!
86  !! mesh
87  call mmg2d_passmesh(mmgmesh, mmgsol, m)
88  !!
89  !! level values at the mesh vertexes
90  call mmg2d_set_solsize(mmgmesh, mmgsol, mmg5_vertex, m%nbNd, mmg5_scalar, ier)
91  do ii=1, m%nbNd
92  call mmg2d_set_scalarsol(mmgsol, val(ii), ii, ier)
93  if (ier/=1) exit
94  end do
95  !!
96  !! clear the input mesh
97  call clear(m)
98 
99  !! #####################
100  !!
101  !! SET MMG PARAMETERS
102  !!
103  !! verbosity, default = 0
104  !!
105  if ( present(verb) ) then
106  call mmg2d_set_iparameter(mmgmesh,mmgsol,&
107  & mmg2d_iparam_verbose, verb, ier)
108  else
109  call mmg2d_set_iparameter(mmgmesh,mmgsol,&
110  & mmg2d_iparam_verbose, 0, ier)
111 
112  end if
113  !!
114  !! Hausdorff distance, default = 0.01
115  !!
116  if ( present(hausd) ) then
117  call mmg2d_set_dparameter(mmgmesh, mmgsol, &
118  & mmg2d_dparam_hausd, hausd, ier)
119  else
120  call mmg2d_set_dparameter(mmgmesh, mmgsol, &
121  & mmg2d_dparam_hausd, 0.01_dp, ier)
122  end if
123  !!
124  !! Hausdorff distance, default = 0.01
125  !!
126  if ( present(hmax) ) then
127  call mmg2d_set_dparameter(mmgmesh, mmgsol, &
128  & mmg2d_dparam_hmax, hmax, ier)
129  end if
130  !!
131  !! ISOVAL
132  !!
133  call mmg2d_set_iparameter(mmgmesh,mmgsol,mmg2d_iparam_iso,1,ier)
134 
135  !! #####################
136  !!
137  !! REMESH
138  !!
139  call mmg2d_mmg2dls(mmgmesh, mmgsol, ier)
140 
141  ! call mmg2d_saveMesh(mmgMesh,"fin.mesh",8, ier)
142  ! call mmg2d_saveSol(mmgMesh,mmgSol,"fin.sol",7, ier)
143 
144  !! #####################
145  !!
146  !! END
147  !!
148  !! retrieve the mesh from MMG
149  call mmg2d_retrievemesh(m, mmgmesh, mmgsol)
150 
151  !! deallocate the mmg structures
152  !!
153  call mmg2d_free_all(mmg5_arg_start, &
154  mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
155  mmg5_arg_end)
156 
157 
158  case(3)
159  if ( m%cell_Count(cell_tet)==0 ) call quit("mesh_mmg: reMesh_levelSet:&
160  & dimension 3 mesh withour tetrahedra")
161 
162  call quit('mesh_mmg: reMesh_levelSet:&
163  & dim 3 to be done')
164 
165  case default
166  call quit("mesh_mmg: reMesh_levelSet: wrong dimension")
167 
168  end select
169 
170 #else
171 
172  !! preproc variable WMMG not declared
173  !!
174  call quit("mesh_mmg: reMesh_levelSet: set 'WITH_MMG=1' in the&
175  & file 'CONFIG.in', re-configure (mmg must be installed)")
176 #endif
177 
178  end subroutine remesh_levelset
179 
180  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
181  !!
182  !!
183  !! HERE STARTS THE CODE INCLUDING MMG DATA / FUNCTIONS
184  !!
185  !!
186  !!
187  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188 #ifdef WMMG
189 
190  subroutine mmg2d_passmesh(mmgMesh, mmgSol, m1)
191  mmg5_data_ptr_t, intent(out) :: mmgmesh, mmgsol
192  type(mesh) , intent(in) :: m1
193 
194  integer :: ier, nNd, nEd, nTr, cl, ii, j1
195 
196  if (choral_verb>0) write(*,*) &
197  & "mesh_mmg : mmg2D_passMesh"
198 
199  mmgmesh = 0
200  mmgsol = 0
201 
202  call mmg2d_init_mesh(mmg5_arg_start, &
203  mmg5_arg_ppmesh, mmgmesh, mmg5_arg_ppmet, mmgsol, &
204  mmg5_arg_end)
205 
206  nnd = m1%nbNd
207  ned = m1%cell_Count(cell_edg)
208  ntr = m1%cell_Count(cell_trg)
209 
210  call mmg2d_set_meshsize(mmgmesh, nnd, ntr, ned, ier)
211 
212  !! define mesh vertexes
213  !!
214  do ii=1, nnd
215  call mmg2d_set_vertex(mmgmesh, &
216  &m1%nd(1,ii), m1%nd(2,ii), 0, ii, ier)
217  end do
218 
219  !! vertex cells --> required vertexes
220  !!
221  do cl=1, m1%nbCl
222  if (m1%clType(cl) == cell_vtx) then
223 
224  !! ii = index of the vertex in m%nd
225  j1 = m1%clToNd%row(cl)
226  ii = m1%clToNd%col(j1)
227 
228  call mmg2d_set_vertex(mmgmesh, &
229  &m1%nd(1,ii), m1%nd(2,ii), m1%clTag(cl), ii, ier)
230 
231  call mmg2d_set_requiredvertex(mmgmesh, ii, ier)
232 
233  end if
234  end do
235 
236 
237  !! define mesh edges
238  !!
239  ii = 0
240  do cl=1, m1%nbCl
241  if (m1%clType(cl) == cell_edg) then
242  ii = ii + 1
243  j1 = m1%clToNd%row(cl)
244 
245  call mmg2d_set_edge(mmgmesh, &
246  & m1%clToNd%col(j1), m1%clToNd%col(j1+1), &
247  & m1%clTag(cl), ii, ier)
248 
249  end if
250  end do
251 
252  !! define mesh triangles
253  ii = 0
254  do cl=1, m1%nbCl
255  if (m1%clType(cl) == cell_trg) then
256  ii = ii + 1
257  j1 = m1%clToNd%row(cl)
258 
259  call mmg2d_set_triangle(mmgmesh, &
260  & m1%clToNd%col(j1), m1%clToNd%col(j1+1), m1%clToNd%col(j1+2),&
261  & m1%clTag(cl), ii, ier)
262  end if
263  end do
264 
265  end subroutine mmg2d_passmesh
266 
267  subroutine mmg2d_retrievemesh(m2, mmgMesh, mmgSol)
268  type(mesh) , intent(inout) :: m2
269  mmg5_data_ptr_t, intent(inOut) :: mmgmesh, mmgsol
270 
271  logical, dimension(:) , allocatable :: vtx_req
272  integer, dimension(:) , allocatable :: vtx_ref
273  integer, dimension(:,:), allocatable :: clTab
274 
275  integer :: ier, nNd, nEd, nTr, nb_vtx_cell, ii, jj, ll
276  integer :: v1, v2, v3, req, ref
277  real(DP) :: x, y
278 
279  if (choral_verb>0) write(*,*) &
280  & "mesh_mmg : mmg2D_retrieveMesh"
281 
282  call clear(m2)
283 
284  call mmg2d_get_meshsize(mmgmesh, nnd, ntr, ned, ier)
285 
286  call allocmem(m2%nd, 3, nnd)
287  call allocmem(vtx_req, nnd)
288  call allocmem(vtx_ref, nnd)
289  nb_vtx_cell = 0
290  do ii=1, nnd
291  call mmg2d_get_vertex(mmgmesh, x, y, vtx_ref(ii), ll, req, ier)
292  m2%nd(1,ii) = x
293  m2%nd(2,ii) = y
294  m2%nd(3,ii) = 0.0_rp
295 
296  if (req/=0) then
297  vtx_req(ii) = .true.
298  nb_vtx_cell = nb_vtx_cell + 1
299  else
300  vtx_req(ii) = .false.
301  end if
302  end do
303 
304  !! total number of cells = required vertexes + edges + triangles
305  ii = nb_vtx_cell + ned + ntr
306 
307  !! allocations to build the mesh m2
308  call allocmem(m2%clType, ii)
309  call allocmem(m2%clTag , ii)
310  call allocmem(cltab, 3, ii )
311 
312  !! set cell types
313  ii = nb_vtx_cell
314  m2%clType(1:ii) = cell_vtx
315  jj = ii + ned
316  m2%clType(ii+1:jj) = cell_edg
317  m2%clType(jj+1:) = cell_trg
318 
319  !! set cell nodes
320  !!
321  !! 1- vertex cells
322  !!
323  jj = 0
324  do ii=1, nnd
325  if (vtx_req(ii)) then
326  jj = jj + 1
327  cltab(1, jj) = ii
328  m2%clTag(jj) = vtx_ref(ii)
329  end if
330  end do
331  call freemem(vtx_req)
332  call freemem(vtx_ref)
333  !!
334  !! 2- edge cells
335  !!
336  jj = nb_vtx_cell
337  do ii=1, ned
338 
339  !! index of the edge ii in the mesh m2
340  jj = jj + 1
341 
342  call mmg2d_get_edge(mmgmesh, v1, v2, ref, ll, req, ier)
343  cltab(1, jj) = v1
344  cltab(2, jj) = v2
345  m2%clTag(jj) = ref
346 
347  end do
348  !!
349  !! 3- triangle cells
350  !!
351  jj = nb_vtx_cell + ned
352  do ii=1, ntr
353 
354  !! index of the edge ii in the mesh m2
355  jj = jj + 1
356 
357  call mmg2d_get_triangle(mmgmesh, v1, v2, v3, ref, req, ier)
358  cltab(1, jj) = v1
359  cltab(2, jj) = v2
360  cltab(3, jj) = v3
361  m2%clTag(jj) = ref
362  end do
363 
364  !! end the mesh m2 construction
365  !!
366  call mesh_create_cltab(m2, cltab)
367  call mesh_create_end(m2)
368 
369  if (.NOT.valid(m2)) call quit('mesh_mmg: &
370  & mmg2d_retrieveMesh: construction not valid')
371 
372  end subroutine mmg2d_retrievemesh
373 
374 #endif
375 
376 end module mesh_mmg
Derived type mesh.
Definition: mesh_mod.F90:92
is the csr matrix well defined ?
Definition: csr_mod.F90:101
DERIVED TYPE graph: sparse matrices of boolean in CSR format
Definition: graph_mod.F90:33
HANDLE THE COMMUNICATION WITH THE LIBRARY MMG
Definition: mesh_mmg.F90:17
integer, parameter cell_tet
Tetrahedron.
BASIC TOOLS
Definition: basic_tools.F90:8
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
subroutine, public quit(message)
Stop program execution, display an error messahe.
integer, parameter cell_edg
Edge (line segment)
destructor
Definition: real_type.F90:127
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
integer, parameter cell_trg
Triangle.
CHORAL CONSTANTS
DERIVED TYPE mesh
Definition: mesh_mod.F90:57
allocate memory for real(RP) arrays
Definition: real_type.F90:158
integer choral_verb
Verbosity level.
subroutine, public mesh_create_end(m)
mesh constructor : final step
Definition: mesh_mod.F90:567
subroutine, public mesh_create_cltab(m, clTab)
mesh constructor, second step, cell types and cell nodes are known cell types are in mclType cell nod...
Definition: mesh_mod.F90:380
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
integer, parameter cell_vtx
Vertex.
DEFINITION OF GEOMETRICAL CELLS (for meshes)
Definition: cell_mod.f90:82
subroutine, public remesh_levelset(m, val, verb, hausd, hMax)
Definition: mesh_mmg.F90:54