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