swh:1:snp:6d9b0a7121ef7aee4569095d3304fd3c77f27314
Raw File
Tip revision: 04a7cf3445bc390bc8dfaf93341a1f8f2008072c authored by Software Heritage on 10 April 2020, 08:52:33 UTC
hal: Deposit 544 in collection hal
Tip revision: 04a7cf3
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
back to top