Poissons Equation with Point source

Hello everyone :)

in matlab´s Example "Poisson's Equation with Point Source and Adaptive Mesh Refinement" ist a function "circlef" mentioned, that create the point source by returning 1/area for the triangle that contains the origin and zero elsewhere. How do I create such a funtion? Meaning how can I define f differently on specific traingles in the mesh?

Thanks in advance

Hannah

[Dyuman Joshi]

"How do I create such a funtion?"

Use the description given as its definition.

"Meaning how can I define f differently on specific traingles in the mesh?"

Use if, elseif, else conditions.

[Hannah Burmester]

Thanks you for your reply. Yes I thought of something like

if (0,0) in triangle

f = 1/area(triangle)

else

f = 0

end

But how can I translate that into code?

Answer 1

[Dyuman Joshi]

You can use inpolygon to see if the origin is inside the triangle or not and polyarea to calculate the area of the triangle -

Note - inpolygon() also considers the points that lie ON the edges of the polygon (see the case below).

x = [-2 -1 1];
y = [-2 1 -1];
plot(polyshape(x,y))
hold on
plot(0, 0, '.r', 'MarkerSize', 10)

if inpolygon(0, 0, x, y)
    f = 1./polyarea(x,y);
else 
    f = 0;
end
f

f = 0.2500

This can be condensed to -

F = inpolygon(0, 0, x, y)./polyarea(x,y)

F = 0.2500

[Dyuman Joshi]

If the criteria you use does not involve the points ON the edge to be counted, use this -

x = [-2 -1 1];
y = [-2 1 -1];


[in,on] = inpolygon(0, 0, x, y)

in = logical
   1

on = logical
   1

if in & ~on
    f = 1./polyarea(x,y);
else 
    f = 0;
end
f

f = 0

Which, again, can be condensed to -

F = (in & ~on)./polyarea(x,y)

F = 0

[Dyuman Joshi]

Any updates, @Hannah Burmester?

[Torsten]

I'd simply edit the mentionned MATLAB function:

edit circlef

Isn't that possible ?

[Hannah Burmester]

Thank you all for your replays! Editing the matlab function should work.

[Dyuman Joshi]

From a cursory glance, I'd say my solution is simpler than whatever is being done here -

type circlef.m

function f=circlef(p,t,u,time)

%point source at (0,0)

% Copyright 1994-2001 The MathWorks, Inc.

x=0;
y=0;

np=size(p,2);
nt=size(t,2);

[ar,t1,t2,t3]=pdetrg(p,t);

[t1,tn,t2,t3]=tri2grid(p,t,zeros(np,1),x,y);

f=zeros(1,nt);

if ~isnan(tn)
  f(tn)=1/ar(tn);
end