code does not call a designed function that works in other program

Hello and sorry for this long query. I am not very good at this , I have to admit it.

I am encountering an issue when trying to apply custom boundary conditions in a PDE model using MATLAB. Specifically, I have a custom boundary condition function mybc1 that works correctly in another code when using structuralBC to apply displacement conditions to vertices. However, in my current code, when I try to use applyBoundaryCondition to set mixed boundary conditions on specific faces, it seems that the solver does not enter the custom boundary condition function mybc1.

Here's a brief description of what I am trying to do:

1. Model Definition: I have defined a PDE model with certain geometry and mesh.
2. Boundary Condition Application: I attempt to apply mixed boundary conditions on specific faces using applyBoundaryCondition with the u parameter set to a custom function handle @(location, state)mybc1(location, state, dis, ts, in).
3. Issue: It appears that the solver does not call the mybc1 function during the solution process. I verified this by setting breakpoints inside mybc1, which are never hit.

The mybc1 function interpolates the acceleration data (dis) at the times specified by the solver (state.time) and returns the corresponding values. This function has been tested and works correctly in another program with structuralBC.

code as follows:

%____________Define dimensions of the plate________________________________
len = 1.22;                          % length in x direction
width = 1.22;                        % Width in y direction
sq_side = 0.10;                      % Side length of the squares in corners
% Create the geometry description matrix___________________________________
% -----Define the outer boundary-------------------------------------------
outer_boundary = [3, 4, 0, len, len, 0, 0, 0, width, width]';
%------Define squares in the corners---------------------------------------
square_1 = [3, 4, 0, sq_side, sq_side, 0, 0, 0, sq_side, sq_side]';
square_2 = [3, 4, len - sq_side, len, len, len - sq_side, 0, 0, sq_side, sq_side]';
square_3 = [3, 4, len - sq_side, len, len, len - sq_side, width - sq_side, width - sq_side, width, width]';
square_4 = [3, 4, 0, sq_side, sq_side, 0, width - sq_side, width - sq_side, width, width]';
%------Combine all the geometries------------------------------------------
gdm = [outer_boundary, square_1, square_2, square_3, square_4];
%------Define the names for each region------------------------------------
ns = (char('R1','R2','R3','R4','R5'))';
sf = 'R1+R2+R3+R4+R5';
%-------Create the geometry------------------------------------------------
model = createpde(2);
geometryFromEdges(model, decsg(gdm, sf, ns));
%-------Generate the mesh and plot the geometry----------------------------
msh = generateMesh(model, 'Hmax', 0.1);  % Use 'Hmax' to control mesh density
figure
pdemesh(model)

%_____________Create the different regions of the model___________________
% Obtain nodes and elements of the mesh
[p,e,t] = meshToPet(model.Mesh);
% ----------Plot with Faces labels----------------------------------------
figure;
pdegplot(model, 'FaceLabels', 'on', 'FaceAlpha', 0.5);

figure
pdemesh(model, 'NodeLabels', 'on');

% Definition of the constants
E = 4E9;
h_thick = 0.05;
nu = 0.3;
mass = 100;
D = E*h_thick^3/(12*(1-nu)^2);
% 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;
%------Display the symbolic coefficients-----------------------------------
structfun(@disp,symcoeffs);

symcoeffs=subs(symcoeffs,pres,1);
coeffs=pdeCoefficientsToDouble(symcoeffs);
%------pass these coefficients to the pde model----------------------------
specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d,'c',...
    coeffs.c,'a',coeffs.a,'f',coeffs.f);
% INITIAL CONDITIONS ------------------------------------------------------
setInitialConditions(model,[0;0],[0;0]);
%USE THE DATA FROM ASHMOLEAN AS BOUNDARY CONDITIONS------------------------
%load the vector containing the acceleration values
load(totaldata.mat', 'DATA',...
    'TIME','fs');

Unable to resolve the name 'totaldata.mat'.

ac=DATA;
% -------Create vector time, not starting in 0
num_samples = size(DATA, 1);
dt = 1 / fs;                            % interval between measurements
t = (1:num_samples) * dt;               % vector of times
% TRANSFORM ACCELERATIONS INTO DISPLACEMENTS_______________________________
acvedi=AccVelDis(962.54,1005.11,992.35,ac,fs,t);
dis=acvedi{3};
t=[0,t];                                  % I start the time at 0
% Boundary conditions in the faces using the displacement vectors
in=0;                                                                      %counter variable
Faces=[1,2,3,5];                                                           %Vector with faces numbers
for jk = 1:numel(Faces)
    in=in+1;%add 1 in every loop to access the different measurmentes contained in the matrix acceleration
    face= Faces(jk);
    dis1=dis(:,in);
    
    applyBoundaryCondition(model,"mixed","Face",[face,face],"u",@(location,state)mybc1(location,state,dis1,t),"EquationIndex",1,"q",[0 0],"g",0); 
    %applyBoundaryCondition(model,"dirichlet","Face",face,"h",[0 0 ; 0 1],"r",bcfunc)
end
in=0;
tim=[t(2) t(3)];                                                           %tim represents the time of interest to solve the pde
res=solvepde(model, tim);
% Access the solution at the nodal locations
sol=res.NodalSolution;

And this is the function that works in another code:

function bcMatrix = mybc1(location,state,acs,ti)
%This function is use to extrapolate the values of acceleration for the
%time instances chosen by the system to do the integration of the system.
%as imput it requires the measurements matrix and the corresponding time
%measurements. As output, it provices the extrapoleted values for the state
%time

    T=state.time;      
    % Check if T is NaN and assign the previous value if true
    if isnan(T)
       vq = NaN(size(location.x));     %this is what documentation says????     
    else
        ac = acs;                      %gets the value of the corresponding acceleration
        vq = interp1(ti, ac, T);       %interpolates it
    end  
    bcMatrix = vq;
end  

I would really appreciate any help/guidance

[Jorge Garcia Garcia]

DATA is a matrix of 30000x4 and TIME a time vector of 30000x1, fs is the sampling frequency=2048

[Torsten]

I think MATLAB ignores inputs for "u" if you define the type of boundary condition as "mixed". "u" is reserved for Dirichlet boundary conditions.

And why do you repeat the face for which you want to set the condition: [face face] ?

[Walter Roberson]

load(totaldata.mat', 'DATA',...
    'TIME','fs');

bad quoting: you need

load('totaldata.mat', 'DATA',...
    'TIME','fs');

Unfortunately you did not happen to include that file for us to test with.

[Sergio E. Obando]

You have a 2D problem, I think the boundary conditions should be set at the edges rather than the faces. Not all edges in F1-F5 are exterior boundaries.

[Jorge Garcia Garcia]

Hello and thanks for the answer.

I try with edges and seem to be working. Well it works for a limited length of the vector time (if I get the vector with 3-10 values of time).When i try to calculate for a vector with most of time values, I get a mistake:

funny enough if i pass the time as tim= [t(1) t(2) t(3) t(4) t(5)]; it fails, but if I pass time as tim= [t(2) t(6) t(9) t(10) t(15)]; it works

Error using daeic3

Need a better guess y0 for consistent initial conditions.

Error in ode15s (line 301)

[y,yp,f0,dfdy,nFE,nPD,Jfac,dMfac] = daeic3(ode,odeArgs,tspan,htry,Mtype,Mt,Mfun,...

Error in pde.internal.FESolverUtilities/solveTimeDependent (line 102)

sol=ode15s(fcn,tlist,uu0,odeoptions);

Error in pde.PDEModel/solvepde (line 57)

[u,dudt] = self.solveTimeDependent(coefstruct, u0, ut0, tlist, ...

Error in walldividedinfacesinthecorner_31_07_24 (line 116)

res=solvepde(model, tim);

and not really know what to do with it, as boundaries need to be setInitialConditions [0;0], [0;0]

I attached the matrix but it is much shorter that it needs to be due to allowed size requirements. Thanks again

the other function i use is

function [outputs]= AccVelDis(sens0,sens1,sens2,rawaccs,fs,TIME,varargin)
%AccVelDis transforms the raw acceleration in filtered accelerations, velocities and
%displacements. The inputs are the sensititivies of the accelerometers in
%channel 0,1,2,the raw matrix of the acceleration recorded, sampling
%frequency (fs)
%an optional argument is sens3 as sometimes the four chanels are not
%recording. This is past at the end as varagin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------Check if there is any value for sens3-----------------------------                
if length(varargin) >= 1
    sens3 = varargin{1};
else
    sens3 = [];
end
sensit0=sens0/1000;                %Transforms the signal channel 0 mvolts/g into volts/g
sensit1=sens1/1000;                %Transforms the signal channel 1 mvolts/g into volts/g    
sensit2=sens2/1000;
sensit3 = [];
% check that the value is not empty
if ~isempty(sens3)
    sensit3 = sens3 / 1000;
end

[rows, cols] = size(rawaccs);
acceloriginal = zeros(rows, cols);
acceloriginal(:,1)=9.8*rawaccs(:,1)/sensit0;%with this we transform the signal first in g and then in m/s^2
acceloriginal(:,2)=9.9*rawaccs(:,2)/sensit1;%with this we transform the signal first in g and then in m/s^2
acceloriginal(:,3)=9.8*rawaccs(:,3)/sensit2;%with this we transform the signal first in g and then in m/s^2
% Check that there is a 4th column in rawaccs for chan4 and sens3 is not
% empty
if ~isempty(sens3) && cols >= 4 && ~isempty(rawaccs(:, 4))
    % Transform the signal first in g and then in m/s^2
    acceloriginal(:, 4) = 9.8 * rawaccs(:, 4) / sensit3;
end
t=TIME;
%--------------Create vector with 0----------------------------------------
num_zeros = numel(acceloriginal(:,1));
[Numfil,Numcol]=size(acceloriginal);
%zeros_vect = zeros(numel(acceloriginal(:,1)),numel(acceloriginal(1,:)));
zeros_vect = zeros(Numfil, Numcol);
%--------------Pack the matrix with 0 before and after the
%values.---------that

packeddata=[zeros_vect;acceloriginal;zeros_vect];
%------------Filter parameters---------------------------------------------

Fs = fs; % Frecuencia de muestreo en Hz (ajustar según tu configuración)
Fpass = 12;%8;%10; % Frecuencia de paso en Hz (puedes ajustarla según tus necesidades)
Fstop = 9;%5;%5; % Frecuencia de rechazo en Hz
Apass = 1; % Ondulación máxima en la banda de paso (dB)
Astop = 60; % Atenuación mínima en la banda de rechazo (dB)

%--------------Design filter-----------------------------------------------
filterb= designfilt('highpassiir', ...
    'StopbandFrequency', Fstop, ...
    'PassbandFrequency', Fpass, ...
    'StopbandAttenuation', Astop, ...
    'PassbandRipple', Apass, ...
    'SampleRate', Fs);


%filterb=designfilt('highpassiir','StopbandFrequency',10,'PassbandFrequency',10.5,'StopbandAttenuation',60,'PassbandRipple',1,'SampleRate',fs);
%--------------Filter accs first-------------------------------------------
filteraccs =filtfilt(filterb,packeddata);


%-----------VELOCITIES CALCULATIONS---------------------------------------
%---------------Remove the packing of 0?
filteraccs=filteraccs(num_zeros+1:end-num_zeros,:);
vels= cumtrapz(t, filteraccs,1);
%----------Check accelerations with graphs------------------------------
% First subplot
% subplot(3, 1, 1);
% plot(t, acceloriginal(:,1));
% hold on;
% plot(t, filteraccs(:,1));
% title('Raw acc chanel 0 and filter acc');
% ylabel('Acceleration (m/s^2)');
% xlabel('time (s)');

% Second subplot
% subplot(3, 1, 2);
% plot(t, acceloriginal(:,2));
% hold on;
% plot(t, filteraccs(:,2));
% title('Raw acc chanel 1 and filter acc');
% ylabel('Acceleration (m/s^2)');
% xlabel('time (s)');;

% Third subplot
% subplot(3, 1, 3);
% plot(t, acceloriginal(:,3));
% hold on;
% plot(t, filteraccs(:,3));
% title('Raw acc chanel 2 and filter acc');
% ylabel('Acceleration (m/s^2)');
% xlabel('time (s)');
% figure
% periodogram(filteraccs(:,1),[],[],fs);hold on;
% periodogram(filteraccs(:,2),[],[],fs);hold on;
% periodogram(filteraccs(:,3),[],[],fs);


%---------------Filter the velocities------------------------------------
%before filtering I need to pack with 0
packedvels=[zeros_vect;vels;zeros_vect];
filtervels =filtfilt(filterb,packedvels);
%remove the packing of 0
filtervels=filtervels(num_zeros+1:end-num_zeros,:);
%  Integrate velocities to displacements
dis= cumtrapz(t, filtervels,1);
% Pack with 0 to filter------------------------------------------------
packeddis=[zeros_vect;dis;zeros_vect];
% filter dis---------------------------------------------------------------
filterdis =filtfilt(filterb,packeddis);
filterdis=filterdis(num_zeros+1:end-num_zeros,:);
%storing the outputs
% adding a 0 at the start. As 0 t measurements = 0

filteraccsnew=zeros(rows+1,cols);filteraccsnew(2:end, :)=filteraccs;
filtervelsnew=zeros(rows+1,cols);filtervelsnew(2:end, :)=filtervels;
filterdisnew=zeros(rows+1,cols);filterdisnew(2:end, :) = filterdis;
outputs = {filteraccsnew, filtervelsnew, filterdisnew};
end

[Sergio E. Obando]

If there is an issue around a particular time step, you could use a coarse time step for certain parts of the solution. E.g. say the issue is around tk = 5; then run the problem tI_1= 0:0.1:4.9., and then use the solution as the initial conditions and solve with tI_2 = 4.9:0.01:6

See Initial Conditions is Prev Solution. Hope that helps.

[Jorge Garcia Garcia]

hello again!

Back after some holidays and this thing keeps on hating me

The never ending problem with this: Now I am trying to interpolate solution in st points: lets say: 0.5,0.5

I am trying

intervalue=interpolateSolution(res, 0.5, 0.5, :,1); as I want all the time behaviour of that point corresponding to the equation 1, that will give me the displacements. and it says: Unrecognized function or variable 'interpolateSolution'. I am extremelly confused.

[Steven Lord]

The documentation page for the interpolateSolution function states that the first input argument must be "specified as a StationaryResults object, a TimeDependentResults object, or an EigenResults object". Can you confirm that res is one of those types of objects? What does these commands show when run on the res variable?

class(res)
isa(res, 'pde.StationaryResults')
isa(res, 'pde.TimeDependentResults')
isa(res, 'pde.EigenResults')

If none of those isa calls return true, your res variable isn't of a type that interpolateSolution can handle. For example, if you tried to call interpolateSolution(1, 2, 3) you'd likely receive the same type of error.

[Torsten]

intervalue=interpolateSolution(res, 0.5, 0.5, :,1);

I doubt that : as fourth argument is a legal call to interpolateSolution.

[Jorge Garcia Garcia]

Hello and thanks for your answers: I can confirm that res is a TimeDependentResults.

I have tried also with:

-ptoizdo=interpolateSolution(res,0.5,0.5)

Error using pde.PDEResults.interpolateSolutionInternal

Incorrect number of input arguments.

Error in pde.TimeDependentResults/interpolateSolution (line 82)

uintrp = pde.PDEResults.interpolateSolutionInternal(obj,varargin{:});

- ptoizdo=interpolateSolution(res,0.5,0.5,:)

Unrecognized function or variable 'interpolateSolution'.

[Torsten]

Look at the example "Interpolate Time Dependent System" under

https://uk.mathworks.com/help/pde/ug/pde.stationaryresults.interpolatesolution.html#buwvykp-1

on how to call "'interpolateSolution" for a time-dependent system.

Here is the essence:

yy = 0:0.2:10;
zz = 0.5*ones(size(yy));
xx = 10*zz;
component = 3;
uintrp = interpolateSolution(results,xx,yy,zz, ...
                   component,1:length(tlist));

Since your problem is 2d, zz should not appear in your input list to "interpolateSolution".

[Jorge Garcia Garcia]

Thanks!! I manage to access the info:

interpolate(posx,poxy,variable number, timevector)