for loop and function handle

Hello, I would really appreciate an answer as this is driving me crazy.

I have a system of PDE, and I have a plate geometry. I have stored the values of the edges in the perimeter in one vector and the values of the inner edges in another vector.

When I apply boundary conditions to the perimeter:

in=0;   
for k = 1:numel(perimeter)
    edge = perimeter(k);
    in=in+1;
    bcfunc= @(location,state) mybc(location,state,dis,ts,in);
    applyBoundaryCondition(modelTwoDomain,"mixed","Edge",[edge,edge],"u",bcfunc,"EquationIndex",1,"q",[0 0],"g",0); 
end
in=0;

where

function bcMatrix = mybc(location,state,acs,ti,i)

T=state.time; 
% Check if T is NaN and assign the previous value if true
if isnan(T)
    %T = previousTime;
    vq = NaN(size(location.x));    
else
    %previousTime = T;  % Update previousTime with the current value of T
    ac = acs(i,:);
    vq = interp1(ti, ac, T);
    
end

bcMatrix = vq;  
end

it works. I pass dis, a matrix containing rows that represent the displacement of every edge in the perimeter, ts the time measurement and in and index to move along the rows depending on the edge.

Now instead of the perimeter I want to do the same in the inner edges, the values stored in inedges

in=0;   
for k = 1:numel(inedges)
    edge = inedges(k);
    in=in+1;
    bcfunc= @(location,state) mybc(location,state,dis,ts,in);
    applyBoundaryCondition(modelTwoDomain,"mixed","Edge",[edge,edge],"u",bcfunc,"EquationIndex",1,"q",[0 0],"g",0); 
end
in=0;

and in this ocasion, the corresponding edges do not enter into mybc, so it does not apply the boundary conditions. Am I missing something here?This is driving me crazy as it does not make sense

Answer 1

[ProblemSolver]

Hello @Jorge Garcia Garcia --

To my understanding of the code snippet that you provided it can be couple of things. To begin with you can try these things:

• Verify if the 'inedges' vector is correct. You can use this before each loop.

disp(inedges); % Print the contents of the inedges vector

• Print a statement of just show if the loop is correct or not something like this. In this way you which edge is being currntly processed and if the entry is correct or not.

for k = 1:numel(inedges)
    edge = inedges(k);
    disp(edge); % Print the current edge value
    in = in + 1;
    bcfunc = @(location,state) mybc(location,state,dis,ts,in);
    applyBoundaryCondition(modelTwoDomain,"mixed","Edge",[edge,edge],"u",bcfunc,"EquationIndex",1,"q",[0 0],"g",0); 
end

• You can either use MATLAB debugger or this snipped provided below to check if your function is providing the desired answers or not. If not, then you know where you are making the mistake.

function bcMatrix = mybc(location,state,acs,ti,i)
    T = state.time;
    
    % Print the inputs for debugging
    disp(location);
    disp(state);
    disp(acs);
    disp(ti);
    disp(i);
    
    % Rest of your code...
    
    bcMatrix = vq;  
end

Based on the details you provided, I could only help you so much. I hope these suggestions help you to figure out where you are going wrong. You should also read this debugging documentation, so that you can apply them as required:

https://www.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html

[Jorge Garcia Garcia]

I can provide the rest of the code, but most of it is to construct the geometry.

I have checked that the edges number contained in inedge are correct, and it loops correctly, it just does not enter the function mybc, whereas if I use the values contained in perimeter, it does enter mybc function. I will try a couple of things tomorrow, and will keep you posted. Will also follow you advise. Thanks

[Jorge Garcia Garcia]

Update. Instead of using the time steps of the solver I imput a time in the function

if I change the code to:

for k = 1:numel(edges)

edge = edges(k);

dis1=dis(k,:); applyBoundaryCondition(modelTwoDomain,"mixed","Edge",edge,"u",mybc(dis1,ts,ts(20)),"EquationIndex",1);

end

res=solvepde(modelTwoDomain,tim);

sol=res.NodalSolution;

Where mybc=

function bcMatrix = mybc(acs,ti,ts)

T=ts;

% Check if T is NaN and assign the previous value if true

if isnan(T)

%T = previousTime;

vq = NaN(size(location.x)); %this is what documentation says????

else

%previousTime = T; % Update previousTime with the current value of T

ac = acs;

vq = interp1(ti, ac, T);

end

bcMatrix = vq;

end

It access the function mybc. But I am interested in the time steps of the solver to extrapolate the solutions so when I change this to :

for k = 1:numel(edges)

edge = edges(k);

dis1=dis(k,:);

applyBoundaryCondition(modelTwoDomain,"mixed","Edge",edge,"u",@(location,state)mybc(location,state,dis1,ts),"EquationIndex",1);

end

res=solvepde(modelTwoDomain,tim);

% Access the solution at the nodal locations

sol=res.NodalSolution;

where: Where mybc=

function bcMatrix = mybc(acs,ti)

T=state.time;

% Check if T is NaN and assign the previous value if true

if isnan(T)

%T = previousTime;

vq = NaN(size(location.x)); %this is what documentation says????

else

%previousTime = T; % Update previousTime with the current value of T

ac = acs;

vq = interp1(ti, ac, T);

end

bcMatrix = vq; #

I put a breakpoint inside of T=state.time, and the flow does not enter the function.

Any suggestion?

[ProblemSolver]

In you updated code:

You are trying to pass the variables dis1 and ts to the function mybc using the anonymous function. Please read this link to understand the definition of function in MATLAB (https://www.mathworks.com/help/matlab/ref/function.html). You have send three input points 'location', 'state', and 'dis1'. Therefore modify the function accordingly.

Secondly, defining the anonymous function (https://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html). PLease read and understand it. You have to define the inuts arguments correctly. The 'mybc' function expects four input arguments, therefore include them as well.

Please post your code properly and questions and explanation on what you have done. Read this: https://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

[Jorge Garcia Garcia]

sorry, it is that so many trials to see what may be wrong, made me copy the wrong code. here is is the whole code to see whether it helps to understand:

%-----------Define dimensions of the plate-------------------------------

len = 3; % length in x direction

width = 2; % Width in y direction

nodes = 1.22; % Distance between the different horizontal and vertical battens

edges=[];

% Create the geometry description matrix-----------------------------------

region_boundaries_y = [0];

dist_y = 0;

interval = nodes;

while dist_y < width

dist_y = dist_y + interval;

if dist_y > width

dist_y = dist_y - (dist_y - width);

elseif abs(dist_y - width) < 1e-10

dist_y = width;

end

region_boundaries_y = [region_boundaries_y dist_y];

end

% Create divisions in X (horizontal divisions) ----------------------------

region_boundaries_x = [0];

dist_x = 0;

interval = nodes;

while dist_x < len

dist_x = dist_x + interval;

if dist_x > len

dist_x = dist_x - (dist_x - len);

elseif abs(dist_x - len) < 1e-10

dist_x = len;

end

region_boundaries_x = [region_boundaries_x dist_x];

end

num_regions_x = length(region_boundaries_x) - 1;

num_regions_y = length(region_boundaries_y) - 1;

%---------Create the PDE Model with a single domain----------------------

model=createpde;

%-CREATE THE GEOMETRY MATRIX---------------------------------------------

% Coordinates of the boundaries in x and y--------------------------------

coords_x = region_boundaries_x';

coords_y = region_boundaries_y';

sq=0.05;

% Create the geometry description matrix

gdm = [];

for k = 1:numel(coords_x)-1

for j = 1:numel(coords_y)-1

% Check if the coordinates are completely inside the external plate

if coords_x(k) >= 0 && coords_x(k) <= len && coords_y(j) >= 0 && coords_y(j) <= width

% Add the vertical dimension

gdm = [gdm; 2, 5, coords_x(k), coords_x(k+1), coords_x(k+1), coords_x(k+1), coords_x(k), coords_y(j), coords_y(j), coords_y(j)+sq, coords_y(j+1), coords_y(j+1)];

end

end

end

%--------------------------------------------------------------------------

gdm=gdm';

%---------------Create the geometry--------------------------------------

[dl,bt] = decsg(gdm);

for j = 1:numel(region_boundaries_x)

for i=1:numel(region_boundaries_y)

indices = find(dl(2,:) == region_boundaries_x(j) & dl(3,:) == region_boundaries_x(j)&...

dl(4,:) == region_boundaries_y(i) & dl(5,:) == region_boundaries_y(i)+sq);

edges = [edges, indices];

end

end

perimeter = find((dl(2,:)== 0) & (dl(3,:) == 0) |...

(dl(4,:) == 0) & (dl(5,:) == 0)|...;

(dl(4,:) == width) & (dl(5,:) == width)|...

(dl(2,:)== len) & (dl(3,:) == len));

%Allows to select just the left edges of the inner squares

edges = setdiff(edges, perimeter);

geometryFromEdges(model,dl);

% Generate the mesh and plot the geometry---------------------------------

msh = generateMesh(model);

%---------Create the different regions of the model-----------------------

% Obtain nodes and elements of the mesh

[p,e,t] = meshToPet(model.Mesh);

modelTwoDomain = createpde(2);

geometryFromEdges(modelTwoDomain,dl);

%

% Now we create the PDE systems as symbolic equations---------------------

syms pres

syms u1(x,y,t) u2(x,y,t)

pdeeq = [-laplacian(u1,[x y])+u2; D*laplacian(u2,[x y])+ mass*diff(u1,t,t)-pres];

symcoeffs = pdeCoefficients(pdeeq,[u1,u2],'Symbolic',true);

c2=symcoeffs.c;

m2=symcoeffs.m;

symcoeffs=subs(symcoeffs,pres,1);

coeffs=pdeCoefficientsToDouble(symcoeffs);

specifyCoefficients(modelTwoDomain,'m',coeffs.m,'d',coeffs.d,'c',...

coeffs.c,'a',coeffs.a,'f',coeffs.f);

% INITIAL CONDITIONS ---------------------------------------------------

setInitialConditions(modelTwoDomain,[0;0],[0;0]);

dis=12x13000 it is this is the dimension of the matrix containing displacements. Basically I want that the first edge in the loop gets the first row of the matrix, and so on.

tim=[ts(20) ts(30) ts(35)];

for k = 1:numel(edges)

edge = edges(k);

dis1=dis(k,:);

applyBoundaryCondition(modelTwoDomain,"mixed","Edge",edge,"u",@(location,state)mybc1(location,state,dis1,ts),"EquationIndex",1);

end

res=solvepde(modelTwoDomain,tim);

% Access the solution at the nodal locations

sol=res.NodalSolution;

In bold is the part of the code that gives me problems. the vector perimeter contains the values: 1 2 7 8 9 13 14 15 16 17 18 19

if I use

for k = 1:numel(perimeter)

edge = perimeter(k);

dis1=dis(k,:);

applyBoundaryCondition(modelTwoDomain,"mixed","Edge",edge,"u",@(location,state)mybc1(location,state,dis1,ts),"EquationIndex",1);

end

where

function bcMatrix = mybc1(location,state,acs,ti)

T=state.time;

if isnan(T)

%T = previousTime;

vq = NaN(size(location.x));

else

ac = acs;

vq = interp1(ti, ac, T);

end

bcMatrix = vq;

end

The code runs without not problem. U stablish a breaking point in the function in space.time, and the code runs into it.

On the other hand if I use what I need:

for k = 1:numel(edges)

edge = edges(k);

dis1=dis(k,:);

applyBoundaryCondition(modelTwoDomain,"mixed","Edge",edge,"u",@(location,state)mybc1(location,state,dis1,ts),"EquationIndex",1);

end

then the code does not enter the function. Vector edges contains the edges I need edges = 3 5 20 22

I really cannot see why in one case it is entering the function mybc1 and in the other case of edges, it is not.