rescale_mesh.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of NS2DDV. %
% %
% Copyright(C) 2011-2018 C. Calgaro (caterina.calgaro@math.univ-lille1.fr) %
% E. Creusé (emmanuel.creuse@math.univ-lille1.fr) %
% T. Goudon (thierry.goudon@inria.fr) %
% A. Mouton (alexandre.mouton@math.univ-lille1.fr) %
% %
% NS2DDV is free software: you can redistribute it and/or modify it under the terms %
% of the GNU General Public License as published by the Free Software Foundation, %
% either version 3 of the License, or (at your option) any later version. %
% %
% NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY %
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A %
% PARTICULAR PURPOSE. See the GNU General Public License for more details. %
% %
% You should have received a copy of the GNU General Public License along with %
% NS2DDV. If not, see <http://www.gnu.org/licenses/>. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[new_p1, new_hmax, new_hmin, NEW_DOMAIN] = rescale_mesh(MESH, OLD_DOMAIN, PARAMETERS)
% INPUT
% MESH The mesh to be rescaled
% OLD_DOMAIN The domain metadata before space rescaling
% PARAMETERS List of simulation parameters
%
% OUTPUT
% new_p1 P1 nodes after space rescaling
% new_hmax Maximal edge length after space rescaling
% new_hmin Minimal edge length after space rescaling
% NEW_DOMAIN The domain metadata after space rescaling
if ~strcmp(OLD_DOMAIN.GEOMETRY, PARAMETERS.DOMAIN.GEOMETRY)
error('Problem while rescaling the mesh: incoherence between old and new geometries');
end
NEW_DOMAIN.GEOMETRY = PARAMETERS.DOMAIN.GEOMETRY;
np = size(MESH.p1,2);
switch PARAMETERS.DOMAIN.GEOMETRY
case {'RECTANGLE', 'DIHEDRON', 'STEP'}
% Rescale the mesh to the unit square
new_p1 = zeros(size(MESH.p1));
new_p1(1,:) = (MESH.p1(1,:)-ones(1,np)*OLD_DOMAIN.XMIN)/(OLD_DOMAIN.XMAX-OLD_DOMAIN.XMIN);
new_p1(2,:) = (MESH.p1(2,:)-ones(1,np)*OLD_DOMAIN.YMIN)/(OLD_DOMAIN.YMAX-OLD_DOMAIN.YMIN);
% Rescale the mesh to the rectangle we want
new_p1(1,:) = new_p1(1,:)*(PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN) ...
+ ones(1,np)*PARAMETERS.DOMAIN.XMIN;
new_p1(2,:) = new_p1(2,:)*(PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YMIN) ...
+ ones(1,np)*PARAMETERS.DOMAIN.YMIN;
if strcmp(PARAMETERS.DOMAIN.GEOMETRY, 'DIHEDRON')
% Rescale the coordinates of dihedral angles
NEW_DOMAIN.XAXIS_ANGLE1 = (OLD_DOMAIN.XAXIS_ANGLE1-OLD_DOMAIN.XMIN) / ...
(OLD_DOMAIN.XMAX-OLD_DOMAIN.XMIN);
NEW_DOMAIN.XAXIS_ANGLE1 = PARAMETERS.DOMAIN.XMIN + ...
NEW_DOMAIN.XAXIS_ANGLE1*(PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN);
NEW_DOMAIN.YAXIS_ANGLE1 = (OLD_DOMAIN.YAXIS_ANGLE1-OLD_DOMAIN.YMIN) / ...
(OLD_DOMAIN.YMAX-OLD_DOMAIN.YMIN);
NEW_DOMAIN.YAXIS_ANGLE1 = PARAMETERS.DOMAIN.YMIN + ...
NEW_DOMAIN.YAXIS_ANGLE1*(PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YMIN);
NEW_DOMAIN.XAXIS_ANGLE2 = (OLD_DOMAIN.XAXIS_ANGLE2-OLD_DOMAIN.XMIN) / ...
(OLD_DOMAIN.XMAX-OLD_DOMAIN.XMIN);
NEW_DOMAIN.XAXIS_ANGLE2 = PARAMETERS.DOMAIN.XMIN + ...
NEW_DOMAIN.XAXIS_ANGLE2*(PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN);
end
if strcmp(PARAMETERS.DOMAIN.GEOMETRY, 'STEP')
NEW_DOMAIN.XSTEP = (OLD_DOMAIN.XSTEP-OLD_DOMAIN.XMIN) / ...
(OLD_DOMAIN.XMAX-OLD_DOMAIN.XMIN);
NEW_DOMAIN.XSTEP = PARAMETERS.DOMAIN.XMIN + ...
NEW_DOMAIN.XSTEP*(PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN);
NEW_DOMAIN.YSTEP = (OLD_DOMAIN.YSTEP-OLD_DOMAIN.YMIN) / ...
(OLD_DOMAIN.YMAX-OLD_DOMAIN.YMIN);
NEW_DOMAIN.YSTEP = PARAMETERS.DOMAIN.YMIN + ...
NEW_DOMAIN.YSTEP*(PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YMIN);
end
% Rescale the bounding box
NEW_DOMAIN.XMIN = PARAMETERS.DOMAIN.XMIN;
NEW_DOMAIN.XMAX = PARAMETERS.DOMAIN.XMAX;
NEW_DOMAIN.YMIN = PARAMETERS.DOMAIN.YMIN;
NEW_DOMAIN.YMAX = PARAMETERS.DOMAIN.YMAX;
case {'CIRCLE'}
% Rescale the mesh to the unit disk
new_p1 = zeros(size(MESH.p1));
new_p1(1,:) = (MESH.p1(1,:)-ones(1,np)*OLD_DOMAIN.XCENTER)/OLD_DOMAIN.RADIUS;
new_p1(2,:) = (MESH.p1(2,:)-ones(1,np)*OLD_DOMAIN.YCENTER)/OLD_DOMAIN.RADIUS;
% Rescale the mesh to the disk we want
new_p1(1,:) = new_p1(1,:)*PARAMETERS.DOMAIN.RADIUS + ones(1,np)*PARAMETERS.DOMAIN.XCENTER;
new_p1(2,:) = new_p1(2,:)*PARAMETERS.DOMAIN.RADIUS + ones(1,np)*PARAMETERS.DOMAIN.YCENTER;
% Rescale the bounding disk
NEW_DOMAIN.XCENTER = PARAMETERS.DOMAIN.XCENTER;
NEW_DOMAIN.YCENTER = PARAMETERS.DOMAIN.YCENTER;
NEW_DOMAIN.RADIUS = PARAMETERS.DOMAIN.RADIUS;
case {'FROM_FILE'}
error('Only pre-built domain geometries can be rescaled (RECTANGLE, CIRCLE, DIHEDRON, STEP)');
otherwise
error('Unknown geometry.');
end
% Compute hmax and hmin as the maximal and minimal length of an edge in the mesh
new_hmax = 0.;
new_hmin = max(max(new_p1))-min(min(new_p1));
for i=1:size(MESH.t1,2)
v1 = new_p1(:,MESH.t1(1,i));
v2 = new_p1(:,MESH.t1(2,i));
v3 = new_p1(:,MESH.t1(3,i));
le1 = sqrt((v1(1)-v2(1))^2+(v1(2)-v2(2))^2);
le2 = sqrt((v2(1)-v3(1))^2+(v2(2)-v3(2))^2);
le3 = sqrt((v3(1)-v1(1))^2+(v3(2)-v1(2))^2);
new_hmax = max([new_hmax,le1,le2,le3]);
new_hmin = min([new_hmin,le1,le2,le3]);
end
end