3D PDE modelling different c coefficients for different spatial regions

조회 수: 1 (최근 30일)
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!

채택된 답변

Ravi Kumar
Ravi Kumar 2020년 4월 7일
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

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Geometry and Mesh에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by