1##  Copyright (C) 2014  Daniel Kraft <d@domob.eu>
2##  GNU Octave level-set package.
3##
4##  This program is free software: you can redistribute it and/or modify
5##  it under the terms of the GNU General Public License as published by
6##  the Free Software Foundation, either version 3 of the License, or
7##  (at your option) any later version.
8##
9##  This program is distributed in the hope that it will be useful,
10##  but WITHOUT ANY WARRANTY; without even the implied warranty of
11##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12##  GNU General Public License for more details.
13##
14##  You should have received a copy of the GNU General Public License
15##  along with this program.  If not, see <http://www.gnu.org/licenses/>.
16
17## -*- texinfo -*-
18## @deftypefn  {Function File} {[@var{mesh}, @var{vmap}] =} ls_build_mesh (@var{geom}, @var{phi}, @var{glob} = false)
19##
20## Build a triangle mesh for the level-set geometry described by
21## @var{geom} and by the level-set function @var{phi}.  Note that
22## @var{geom} must contain absolute coordinates as set by
23## @code{ls_absolute_geom}.
24##
25## If @var{glob} is set to @code{true}, then the mesh will contain
26## the whole hold-all domain and not just the interior of the level-set
27## geometry.  It will, however, be constructed such that the boundary
28## of the level-set domain is resolved precisely, and such that the
29## level-set domain can be described as a subset of the mesh triangles.
30##
31## The mesh points will be made up @emph{both} of grid nodes as well as
32## intersection points.  While all intersection points will appear
33## as mesh nodes, unless @var{glob} is set, only @strong{inner} grid nodes
34## will correspond to mesh points.  The points on the mesh are numbered
35## differently to the points in @var{geom}, although the indices can be
36## translated both ways using @var{vmap} (see below).
37##
38## On the other hand, the boundary edges of the mesh match precisely the
39## edges given in @code{@var{geom}.bedges}.  They have the same indices
40## in both data structures.
41##
42## The returned struct @var{mesh} follows the format of the @code{msh}
43## package, which means that it contains these fields:
44##
45## @table @code
46## @item p
47## 2 x @var{np} matrix which contains the x- and y-coordinates
48## of the mesh points in its two rows.
49##
50## @item t
51## 4 x @var{nt} matrix describing the mesh triangles.
52## The first three rows contain the indices of mesh points making up the
53## triangle, and the fourth row contains the index of the geometrical
54## surface this triangle belongs to.  It is always set to @code{1} for now.
55## The points are ordered @strong{counter-clockwise} (in the coordinate-system
56## interpretation of the grid).
57##
58## @item e
59## 7 x @var{ne} matrix describing the edges of the mesh.  They
60## correspond precisely to @code{@var{geom}.bedges}.
61## The rows of this matrix are:
62##
63## @table @asis
64## @item 1--2
65## Start- and end-point index.  They are ordered such that
66## the level-set domain is @strong{on the left} of the edge.
67##
68## @item 3--4
69## Always @code{0} for compatibility.
70##
71## @item 5
72## The gamma component index (as per @code{@var{geom}.gamma}) that
73## contains this edge.  Will be in the range 1--@code{@var{geom}.gamma.n}.
74##
75## @item 6--7
76## Geometrical surface index of the domain parts to the right
77## and left.  They are set to @code{0} and @code{1} for now.
78## @end table
79## @end table
80##
81## The second output, @var{vmap}, describes the mappings between geometrical
82## node / intersection-point indices in @var{geom} and the indices of the
83## mesh points used in @var{mesh}.  It is a struct with these fields:
84##
85## @table @code
86## @item grid
87## Maps the index of a node on the grid to the index of the
88## corresponding mesh point.  The value is @code{NA} if the mesh
89## does not contain the point.
90##
91## @item ispt
92## Maps the index of an intersection point (according to
93## @code{@var{geom}.ispts}) to the index of the corresponding mesh point.
94##
95## @item mesh
96## Maps the index of a mesh point back to either the grid or intersection-point
97## index.  A @strong{positive} value means that the mesh point is on the
98## grid and has the respective index, while a @strong{negative} value means
99## that it is an intersection point with index equal to the absolute value
100## of the entry.
101## @end table
102##
103## @seealso{ls_find_geometry, ls_absolute_geom, msh2m_structured_mesh}
104## @end deftypefn
105
106function [mesh, vmap] = ls_build_mesh (geom, phi, glob = false)
107  if (nargin () < 2 || nargin () > 3)
108    print_usage ();
109  endif
110
111  if (!isfield (geom.ispts, "coord") || !isfield (geom, "nodes"))
112    error ("GEOM misses absolute coordinates, use ls_absolute_geom first");
113  endif
114
115  % Find inside array.
116  L = geom.dim(1);
117  M = geom.dim(2);
118  inside = ls_inside (phi(:));
119
120  % Call internal routine.
121  [mesh, vmap] ...
122    = __levelset_internal_mesh (L, M, inside, ...
123                                geom.nodes.coord, geom.ispts.coord, ...
124                                geom.elem.nodelist, geom.elem.index.inner, ...
125                                geom.elem.index.outer, geom.elem.index.bdry, ...
126                                geom.internal.bdryelSegments, glob);
127  % FIXME: Should be 0 for outside triangles (in global mesh)?
128  mesh.t(4, :) = 1; % Geometrical component.
129
130  % Set mesh vertexMap entries to NA where they are -1.  This is not already
131  % done in the C++ code since it uses integers instead of doubles
132  % for the index arrays being built up.
133  vmap.grid(vmap.grid < 0) = NA;
134  assert (all (vmap.ispt > 0));
135  assert (all (isfinite (vmap.mesh)));
136
137  % Set boundary table from geom.bedges.  We only have to convert
138  % geometrical ispt indices to mesh vertex indices.
139  mesh.e = NA (7, geom.bedges.n);
140  mesh.e(1 : 2, :) = vmap.ispt(geom.bedges.ispts');
141  mesh.e([3:4, 6], :) = 0; % Reserved and outside on the right.
142  mesh.e(5, :) = geom.bedges.comp; % Boundary component.
143  mesh.e(7, :) = 1; % Domain on the left.
144endfunction
145
146%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147% Tests.
148
149% Test for error conditions.
150%!error <Invalid call to>
151%!  ls_build_mesh (1)
152%!error <Invalid call to>
153%!  ls_build_mesh (1, 2, 3, 4)
154
155% Test for missing absolute coordinates.
156%!error <GEOM misses absolute coordinates, use ls_absolute_geom first>
157%!  x = linspace (-10, 10, 10);
158%!  h = x(2) - x(1);
159%!  [XX, YY] = meshgrid (x, x);
160%!  phi = ls_genbasic (XX, YY, "sphere", [0, 0], 5);
161%!  phi = ls_normalise (phi, h);
162%!  geom = ls_find_geometry (phi, h);
163%!  ls_build_mesh (geom, phi);
164
165% Helper function to verify that two points are in counter-clockwise
166% direction with respect to a third one.
167%!function verifyDirection (a, b, centre)
168%!  v1 = [a - centre; 0];
169%!  v2 = [b - centre; 0];
170%!  x = cross (v1, v2);
171%!  assert (x(1 : 2), zeros (2, 1));
172%!  assert (x(3) > 0);
173%!endfunction
174
175% A function to verify certain expected properties / postconditions
176% for a given mesh.
177%!function verifyMesh (geom, phi, glob)
178%!  [mesh, vmap] = ls_build_mesh (geom, phi, glob);
179%!  nPts = prod (geom.dim);
180%!
181%!  % Verify ordering of triangle points.
182%!  for i = 1 : size (mesh.t, 2)
183%!    verifyDirection (mesh.p(:, mesh.t(1, i)), mesh.p(:, mesh.t(2, i)), ...
184%!                     mesh.p(:, mesh.t(3, i)));
185%!  endfor
186%!
187%!  % Go through all grid points and verify them with the mesh.
188%!  for i = 1 : nPts
189%!    gridInd = vmap.grid(i);
190%!    if (isna (gridInd))
191%!      assert (!glob);
192%!      continue;
193%!    endif
194%!
195%!    assert (vmap.mesh(gridInd), i);
196%!    assert (mesh.p(:, gridInd), geom.nodes.coord(i, :)');
197%!  endfor
198%!
199%!  % Go through all intersection points and verify them.
200%!  for i = 1 : geom.ispts.n
201%!    gridInd = vmap.ispt(i);
202%!    assert (vmap.mesh(gridInd), -i);
203%!    assert (mesh.p(:, gridInd), geom.ispts.coord(i, :)');
204%!  endfor
205%!
206%!  % Verify mesh.e against geom.bedges.
207%!  for i = 1 : geom.bedges.n
208%!    assert (mesh.e(1 : 2, i), vmap.ispt(geom.bedges.ispts(i, :)'));
209%!    assert (mesh.e(5, i), geom.bedges.comp(i));
210%!  endfor
211%!endfunction
212%
213% Verify a given phi, by first constructing the geometry and mesh
214% both in global and local mode as well as trying both phi and -phi.
215%!function verifyPhi (phi, XX, YY, h)
216%!  for s = [-1, 1]
217%!    cur = phi * s;
218%!    cur = ls_normalise (cur, h);
219%!    geom = ls_find_geometry (cur, h);
220%!    geom = ls_absolute_geom (geom, XX, YY);
221%!
222%!    verifyMesh (geom, cur, false);
223%!    verifyMesh (geom, cur, true);
224%!  endfor
225%!endfunction
226
227% Test against some example phis.  We do not use ls_get_tests, since
228% we only want phis where we have a grid, too.
229%!test
230%!  n = 20;
231%!  x = linspace (-10, 10, n);
232%!  h = x(2) - x(1);
233%!  [XX, YY] = meshgrid (x, x);
234%!
235%!  phi = ls_union (ls_genbasic (XX, YY, "sphere", [-6, -6], 3), ...
236%!                  ls_genbasic (XX, YY, "sphere", [6, 6], 3));
237%!  verifyPhi (phi, XX, YY, h);
238%!
239%!  phi = ls_setdiff (ls_genbasic (XX, YY, "sphere", [0, 0], 8), ...
240%!                    ls_genbasic (XX, YY, "sphere", [0, 0], 4));
241%!  verifyPhi (phi, XX, YY, h);
242
243%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
244% Demo.
245
246% Build a mesh both in global and non-global mode and plot the result.
247%!demo
248%!  x = linspace (-10, 10, 11);
249%!  h = x(2) - x(1);
250%!  [XX, YY] = meshgrid (x, x);
251%!
252%!  phi = ls_union (ls_genbasic (XX, YY, "sphere", [5, 5], 4.5), ...
253%!                  ls_genbasic (XX, YY, "sphere", [-2, -2], 4));
254%!  phi = ls_normalise (phi, h);
255%!  geom = ls_find_geometry (phi, h);
256%!  geom = ls_absolute_geom (geom, XX, YY);
257%!
258%!  mesh = ls_build_mesh (geom, phi, false);
259%!  meshGlob = ls_build_mesh (geom, phi, true);
260%!
261%!  figure ();
262%!  hold ("on");
263%!  trimesh (meshGlob.t(1 : 3, :)', meshGlob.p(1, :), meshGlob.p(2, :), "r");
264%!  trimesh (mesh.t(1 : 3, :)', mesh.p(1, :), mesh.p(2, :), "g");
265%!  for i = 1 : size (mesh.e, 2)
266%!    plot (mesh.p(1, mesh.e(1 : 2, i)), mesh.p(2, mesh.e(1 : 2, i)), "b");
267%!  endfor
268%!  hold ("off");
269%!  legend ("Global", "Domain", "Boundary");
270%!  axis ("equal");
271%!  axis ("square");
272