How can I delete a circle from my 2D Meshgrid domain?

With the following code, I want to plot the velocity field but the points (x,y) falling within the filled circle have really high values with respect to the normal ones which are outside.

I want to delete (not consider) from my domain, points (x,y) falling inside the circle. This way, I will be able to scale my velocities accordingly.

U0=1; %Uniform flow velocity

r0=1; %radius of the imaginary sphere

k=0.5*U0*r0^3 %strength of the dipole

c=k*4*pi;

x = -2.05:.1:2.05;

y = -2.05:.1:2.05;

[X,Y] = meshgrid(x,y);

t=atan2(Y,X);

r=sqrt(X.^2+Y.^2);

%Stream function

psi = 0.5*U0*Y.^2.*(1-r0^3./(X.^2+Y.^2).^1.5);

%Velocity potential

phi = U0*X.*(1+0.5*r0^3./((X.^2+Y.^2).^1.5));

% Velocity in polar coordinates

vr = U0.*cos(t).*(1-r0^3./(r.^3));

vt = -U0.*sin(t).*(1+.5*r0^3./(r.^3));

% Velocity in cartesian coordinates

vx = U0 + c./(4*pi*(X.^2 + Y.^2).^(3/2)) - 3*c*X.^2./(4*pi*(X.^2 + Y.^2).^(5/2));

vy = -(3*c*X.*Y)./(4*pi*(X.^2 + Y.^2).^(5/2));

%Pressure coefficient on the sphere surface

theta=0:pi/100:2*pi;

Cp=1-9/4.*sin(theta).^2;

%Plots streamlines

contour(X,Y,psi,20,'k')

xlabel('x'); ylabel('y')

hold on

grid on

line(xlim, [0,0], 'Color', 'r', 'LineWidth', 1); % Draws a line for X axis

line([0,0], ylim,'Color', 'r', 'LineWidth', 1); % Draws a line for Y axis

axis equal

title("Potential flow around a sphere")

%Plots equipotential lines

contour(X,Y,phi,400,'k')

%Plots velocities

factor=.001

quiver(X,Y,vx*factor,vy*factor,'g','AutoScale','off');

filling(r0,'w');

[Joseph Vannel Fotso Fogang]

It plots the following without the manually added filled circle. All the points inside the circle are not physically meaningful, they simply need not to be plotted. I don't know how to do so.

Answer 1

[Aditya Patil]

You can use boolean matrix to select unwanted variables, and set them to zero or nan.

For example, in this case, you can find those X,Y pairs which lie inside unit circle as follows,

flag = [abs(sqrt(X.^2 + Y.^2)) < 1] % We are calculating distance from origin

Then use this flag to remove unwanted data from variables to be plotted, as follows

psi(flag) = nan % or zero, as appropriate

Then plot the data as usual.

Answer 2

[David]

The solution from Aditya Patil is a good start, but the resulting excluded area is a square, rather than a circle.

You need to do the same thing, but apply it to the X,Y meshgrid values, rather than the x,y ranges.

Try this:

~~~~~~~~~~~~~~~~~~~~~

U0=1; %Uniform flow velocity

r0=1; %radius of the imaginary sphere

k=0.5*U0*r0^3; %strength of the dipole

c=k*4*pi;

x = -2.05:.1:2.05;

y = -2.05:.1:2.05;

[X,Y] = meshgrid(x,y);

t=atan2(Y,X);

r=sqrt(X.^2+Y.^2);

%Stream function

psi = 0.5*U0*Y.^2.*(1-r0^3./(X.^2+Y.^2).^1.5);

%Velocity potential

phi = U0*X.*(1+0.5*r0^3./((X.^2+Y.^2).^1.5));

% Velocity in polar coordinates

vr = U0.*cos(t).*(1-r0^3./(r.^3));

vt = -U0.*sin(t).*(1+.5*r0^3./(r.^3));

% Velocity in cartesian coordinates

vx = U0 + c./(4*pi*(X.^2 + Y.^2).^(3/2)) - 3*c*X.^2./(4*pi*(X.^2 + Y.^2).^(5/2));

vy = -(3*c*X.*Y)./(4*pi*(X.^2 + Y.^2).^(5/2));

%Pressure coefficient on the sphere surface

theta=0:pi/100:2*pi;

Cp=1-9/4.*sin(theta).^2;

%Plots streamlines

% ---------- Start exclude points within RADIUS -----------

% Save copies of original grid matrices.

% Excluded points are needed for velocity vectors.

X0 = X;

Y0 = Y;

% Exclude plot points inside circle of R = RADIUS

% NOTE: Inserted after calculations, since the excluded points affect

RADIUS = 1;

flag = [abs(sqrt(X.^2 + Y.^2)) < RADIUS];

X(flag) = nan;

Y(flag) = nan;

% Add circle to mark the boundary of excluded points.

viscircles([0 0], RADIUS, 'Color','k', 'LineWidth', 1, 'LineStyle','--');

hold on;

% ---------- Finish exclude points within RADIUS -----------

contour(X,Y,psi,20,'k')

xlabel('x'); ylabel('y')

hold on

grid on

line(xlim, [0,0], 'Color', 'r', 'LineWidth', 1); % Draws a line for X axis

line([0,0], ylim,'Color', 'r', 'LineWidth', 1); % Draws a line for Y axis

axis equal

title("Potential flow around a sphere")

%Plots equipotential lines

contour(X,Y,phi,400,'k')

%Plots velocities

factor=.001;

% NOTE: Using original points for velocity vectors,

% since interesting velocities are in excluded area.

quiver(X0,Y0,vx*factor,vy*factor,'g','AutoScale','off');

filling(r0,'w');

~~~~~~~~~~~~~~~~~~~~~