3D PDE modelling different c coefficients for different spatial regions

1 vue (au cours des 30 derniers jours)
Hi,
I have a geometry imported into MATLAB shown in the picture below.
For reference, the cube describes tissue bed and the cylindrical holes describe blood vessel with dirichlet pressure boundary conditions.
I am trying to specify c coefficients to change the properties in a spherical region within the cube,
with the coefficient inside and outside this spherical region set to different scalar values.
I know you cannot specify subdomains in the same way as in the 2D pde modeller, but I tried to use an "if statement" to spatially set c coefficients without much luck.
function c = c_region_3D(location,state)
%% Input variables for region which we want to alter
% sphere 1
xc_old1=0.1875;
yc_old1=0.1875;
zc_old1= 0.125;
r_old=0.1;
a=1; % counter for all voxels within cube
for i=1:length(location.x)
for j=1:length(location.y)
for k=1:length(location.z)
if ((location.x(i)-xc_old1).^2 + (location.y(j)-yc_old1).^2 + (location.z(k)-zc_old1).^2<=(r_old).^2)==1
c(1,a) = 3.3929e-13;
a=a+1;
else
c(1,a) = 4.4762e-13;
a=a+1;
end
end
end
end
However, I only seem to be able to provide a c matrix with the length of one side of the cube - how do I therefore reference 3D points to specify the c coefficient in this region?
Thank you!

Réponse acceptée

Ravi Kumar
Ravi Kumar le 7 Avr 2020
Hi Barnaby,
If all you want is different values of c for points within sphere and outside sphere, then you can just do this:
function c = c_region_3D(location,state)
xc_old1=0.1875;
yc_old1=0.1875;
zc_old1= 0.125;
r_old=0.1;
c = zeros(size(location.x));
%find IDs of points in sphere or not and use it for logical indexing.
Index_of_PointsInSphere = (location.x-xc_old1).^2 + (location.y-yc_old1).^2 + (location.z-zc_old1).^2<=(r_old).^2;
Index_of_Points_NOT_InSphere = ~Index_of_PointsInSphere;
c(Index_of_PointsInSphere) = 3.3929e-13;
c(Index_of_Points_NOT_InSphere) = 4.4762e-13;
end
Regards,
Ravi

Plus de réponses (0)

Produits

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by