Plot solution of pde toolbox on a line (or submanifold)

I'm using the pde toolbox to solve a certain elliptic equation in 2D.

Solution is fine, although I do need to plot it along a given line, i.e. to cut a planar slice from the 3D mesh representing the solution.

I can't figure out a way that smartly involves the toolbox functions.

Any help appreciated.

Answer 1

[Bill Greene]

Hi,

Here is one way to create such a plot. Assume you have the point matrix created by the PDE Toolbox mesher, p, and a solution vector, u. The function below will create a plot of that solution along a line defined by the x and y locations of the two end points. My example is for a solution on a unit square and I want a plot along the line (0,.5) to (1,.5). I want to include 25 points in the plot. As you can see, the real work is being done by the TriScatteredInterp function from core MATLAB.

plotAlongLine(p, u, [0,.5], [1,.5], 25);

function plotAlongLine(p, u, xy1, xy2, numpts)
x = linspace(xy1(1),xy2(1),numpts);
y = linspace(xy1(2),xy2(2),numpts);
F = TriScatteredInterp(p(1,:)', p(2,:)', u);
uxy = F(x,y);
figure; plot(x, uxy);
end

Bill

[Alessandro]

Perfectly working! thanks!

[Lukasz]

Hi,

I ask for your help, I am a beginner.

I want to get a graph of temperature changes on the edge of a polygon. Between points with coordinates from (0,0) to (1,3).

On the horizontal axis is the distance from point (0,0), and on the vertical axis is the temperature value.

Luk.

[Torsten]

Does this help as an example ?

clear all;
thermalmodel = createpde();
R1= [3,4,0,3,3,0,0,0,10,10]';
R2= [3,4,3,6,6,3,-5,-5,5,5]';
gd= [R1 R2];
sf= 'R1+R2';
ns = char('R1','R2');
ns = ns';
dl = decsg(gd,sf,ns);
geometryFromEdges(thermalmodel,dl);
pdegplot(thermalmodel,"EdgeLabels","on","FaceLabels","on")
xlim([-1 7])
ylim([-6 11])
axis equal

%% Generate and plot mesh
generateMesh(thermalmodel);
figure
pdemesh(thermalmodel)
title("Mesh with Quadratic Triangular Elements")

%% Apply BCs
% Edge 3 is left edge; Edge 5 is right side; edge 8 is middle; Edges 2 & 6 top; Edges 1 & 4 bottom
% Edges 7 & 9 offset vertical edges
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[3],u=100);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[5],u=20);
%% Apply thermal properties in first region
rho1= 1 ;                                  % 
cp1= 1  ;                                   % 
rhocp1= rho1*cp1 ;                 %
k1= 1   ;                                    % W/mK
%% Apply thermal properties in second region
rho2= 10 ;                                % 
cp2= 10  ;                                 % 
rhocp2= rho2*cp2   ;               %
k2= 10    ;                                 % W/mK
%% Define uniform volumetric heat generation rate
Qgen= 0   ;                   % W/m3                      [1e12]
%% Define coefficients for generic Governing Equation to be solved
f= [Qgen];
specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=f, Face=1);
specifyCoefficients(thermalmodel, m=0, d=rhocp2, c=k2, a=0, f=f, Face=2);
coeffs= thermalmodel.EquationCoefficients;
findCoefficients(coeffs, Face=1);
findCoefficients(coeffs, Face=2);
%% Apply initial condition
setInitialConditions(thermalmodel, 20);
%% Define time limits
tlist = 0:1:50;
%% Solve 
thermalresults= solvepde(thermalmodel, tlist);
% Plot results
sol = thermalresults.NodalSolution;
subplot(2,2,1)
pdeplot(thermalmodel,"XYData",sol(:,50), ...
    "Contour","on",... 
    "ColorMap","jet")

plotAlongLine(thermalresults.Mesh.Nodes, thermalresults.NodalSolution(:,50), [0,.5], [6,.5], 50);

function plotAlongLine(p, u, xy1, xy2, numpts)
  x = linspace(xy1(1),xy2(1),numpts);
  y = linspace(xy1(2),xy2(2),numpts);
  F = TriScatteredInterp(p(1,:)', p(2,:)', u);
  uxy = F(x,y);
  figure; plot(x, uxy);
end

[Lukasz]

Thank you for the answer, but unfortunately, I cannot use this solution.

My case:

In PDE Toolbox, I create geometry, mesh, boundary conditions, and simulate heat flow.

Then, using the export command in PDE Toolbox, I export the data "p t e u" to a matrix.

I want to use a script and the "p t e u" data to make a temperature plot on the boundary line between points P1 (0,0) and P2 (1,3).

The figure is above.

Please give me some advice.

[Torsten]

I don't know what "p t e u" data are.

In order to use the function "plotAlongLine", you need a 2xN matrix p of (x,y) coordinates and an 1xN vector u as solution vector in these points.

[Lukasz]

Thank you for your answer.

The video explains what "p e t" is (in 3:33 min), and u, it is the resulting temperature value.

/ https://www.youtube.com/watch?v=Mg2pML84GGY /

All parameters after export from PDE toolbox, are tables in Workspace and available from the command line.

However, I cannot understand the table structure to make a temperature chart on the polygon edge.

The plotAlongLine function does not work correctly.

[Torsten]

It should work exactly as suggested by @Bill Greene :

plotAlongLine(p, u, [0,0], [1,3], 25);

function plotAlongLine(p, u, xy1, xy2, numpts)
  x = linspace(xy1(1),xy2(1),numpts);
  y = linspace(xy1(2),xy2(2),numpts);
  F = TriScatteredInterp(p(1,:)', p(2,:)', u);
  uxy = F(x,y);
  figure; plot(sqrt(x.^2+y.^2), uxy);
end

If not, please include your code.

[Lukasz]

The use of the plotAlongLine function ends with an error.

Details are in the pdf file.

I also tried another solution using the functions pdeInterpolant and evaluate , and later I get the results.

Here is an example function to get the results.

//-----------

t = linspace(0,1);

x = -1 + 0.5*t;

y = 2 - 0.1*t;

uout = evaluate(F,x,y); % Assumes F exists

plot(t,uout)

//-----------

This is very complicated because I have to use parameter t, and there is surely an easier solution.

[Torsten]

Pass the time step to "plotAlongLine":

timestep = 4;
numpts = 25;
pstart = [0 0];
pend = [1 3];
plotAlongLine(p, u, pstart, pend, numpts, timestep);

function plotAlongLine(p, u, xy1, xy2, numpts, timestep)
  x = linspace(xy1(1),xy2(1),numpts);
  y = linspace(xy1(2),xy2(2),numpts);
  F = TriScatteredInterp(p(1,:)', p(2,:)', u(:,timestep));
  uxy = F(x,y);
  figure; plot(sqrt(x.^2+y.^2), uxy);
end

[Lukasz]

The function works correctly in this setup.

This is a very helpful example for me.

Thank you very much for your help!!!

[Torsten]

I think the application of "plotAlongLine" is more flexible if you take the plotting part out of the function:

timestep = 4;
numpts = 25;
pstart = [0 0];
pend = [1 3];
[x,y,uxy] = plotAlongLine(p, u, pstart, pend, numpts, timestep);
plot(sqrt(x.^2+y.^2),uxy)

function [x,y,uxy] = plotAlongLine(p, u, xy1, xy2, numpts, timestep)
  x = linspace(xy1(1),xy2(1),numpts);
  y = linspace(xy1(2),xy2(2),numpts);
  F = TriScatteredInterp(p(1,:)', p(2,:)', u(:,timestep));
  uxy = F(x,y);
end

[Lukasz]

One more question.

How can I export any variables calculated in the function to an external file, especially the data used to create a chart?

[Torsten]

timestep = 4;
numpts = 25;
pstart = [0 0];
pend = [1 3];
[x,y,uxy] = plotAlongLine(p, u, pstart, pend, numpts, timestep);
plot(sqrt(x.^2+y.^2),uxy)
writematrix([(sqrt(x.^2+y.^2)).',uxy.'],'Results.xls')

function [x,y,uxy] = plotAlongLine(p, u, xy1, xy2, numpts, timestep)
  x = linspace(xy1(1),xy2(1),numpts);
  y = linspace(xy1(2),xy2(2),numpts);
  F = TriScatteredInterp(p(1,:)', p(2,:)', u(:,timestep));
  uxy = F(x,y);
end

Answer 2

[Lukasz]

New task and I found a new problem. :(

Problem with determining the temperature on the line with endpoints at points A(0;0) and B(0.88; 5) - the left side of the polygon.

The simulation is about heat transfer.

In the simple example I asked about in previous posts, which is above, everything works.

Details of the problem are in the pdf file.

Please, I ask for your help.

Luk.