Structured mesh in MATLAB PDE toolbox

7 次查看（过去 30 天）
Kaveh Gharibi
Kaveh Gharibi 2017-11-23
评论： Kaveh Gharibi 2017-11-29
I found a way to create an structured quadratic mesh. They way I do it is as follows. For example if I have a rectangle. I divide it to many squares and then mesh each square to two triangular mesh.
if true
numberOfPDE = 3;
model = createpde(numberOfPDE);
global dl
L=2;
h=1;
dl=.05 ;
dh=dl ;
[xq,yq]=meshgrid( [0:dl:L], [-h/2:dh:h/2]);
n=size(xq,1);
m=size(xq,2);
k=1;
R=[];
for i=1:n-1
for j=1:m-1
R(:,k) =[3,4,xq(i,j),xq(min(i+1,n),j),xq(min(i+1,n),min(m,j+1)),xq(i,min(m,j+1)),yq(i,j),yq(min(i+1,n),j),yq(min(i+1,n),min(m,j+1)),yq(i,min(m,j+1))];
k=k+1;
end
end
end
gm = [R];
for i=1:size(R,2)
str{i} = ( sprintf('R%d',i) );
end
sf=[];
for i=1:size(R,2)
sf=char([str{i},'+',sf]);
end
sf=sf(1:end-1);
delete(model.Geometry)
model.Geometry = []
ns = char(str);
ns = ns';
g= decsg(gm,sf,ns);
geometryFromEdges(model,g);
figure
pdegplot(model,'EdgeLabels','on');
axis equal
title 'Geometry with Edge Labels';
This code creates the geometry, The following is the meshing part.
if true
order='quadratic';
% order='linear';
mesh=generateMesh(model,'GeometricOrder', order ,'Hmax' ,.05 , 'MesherVersion','R2013a','Hgrad',1.01,'JiggleIter',100 );
figure(200)
pdemesh(model,'NodeLabels','off')
[p,e,t] = meshToPet(mesh);
axis equal
end
Then PDE is solved as follows.
if true
a = [0 0 0 0 0 0 0 0 0]';
f = [0 0 0]’;
applyBoundaryCondition(model,'mixed',...
'Edge', 761:780 ,...
'u', 0,'EquationIndex', 3 ,'Vectorized','on' );
applyBoundaryCondition(model,'dirichlet',...
'Edge',[1:40,1601:1609],...
'u', @myConstPsi2Rev,'Vectorized','on' );
applyBoundaryCondition(model,'dirichlet',...
'Edge',1610:1611,...
'u', @myConstPsi3Rev,'Vectorized','on' );
applyBoundaryCondition(model,'dirichlet',...
'Edge',[1612:1620,1561:1600],...
'u', @myConstPsi4Rev,'Vectorized','on' );
C=@ccoefPsiDevi;
specifyCoefficients(model,'m',0,'d',0,'c',C,'a',a,'f',f);
end
if true
format long
model.SolverOptions.MinStep=0;
model.SolverOptions.ReportStatistics = 'on';
model.SolverOptions.MaxIterations=5000;
model.SolverOptions.ResidualTolerance=1.6819e-06
%%solve PDE
tic
result = solvepde(model);
toc
end
The functions used as C coefficients and Boundary conditions are as follows.
if true
function cmatrix = ccoefPsiDevi(region,state)
% global xe ye %epsilon %error
n1 = 36;
nr = numel(region.x) ;
cmatrix = zeros(n1,nr);
e= (.1);
normrev= (sqrt(e^2 +(state.uy(1,:)+state.ux(2,:)).^2+2*state.uy(2,:).^2+2*state.ux(1,:).^2 ) );
fac=1 ;
cmatrix(1,:) = fac*2./normrev;
cmatrix(4,:) = fac*1./normrev;
cmatrix(14,:) = fac*1./normrev;
cmatrix(26,:) = -region.y;
cmatrix(27,:) = region.y;
cmatrix(7,:) = fac*1./normrev;
cmatrix(17,:) = fac*1./normrev;
cmatrix(20,:) = fac*2./normrev;
cmatrix(30,:) = region.x;
cmatrix(31,:) = -region.x;
cmatrix(10,:) = -region.y;
cmatrix(11,:) = region.y;
cmatrix(22,:) = region.x;
cmatrix(23,:) = -region.x;
end
end
if true
function bcMatrix = myConstPsi2Rev(region,state)
bcMatrix = [ 0*region.y ;
20*dl*ones(size(region.y));
0*region.y] ;
end
end
if true
function bcMatrix = myConstPsi3Rev(region,state)
global psixb psiyb
bcMatrix = [ 0*region.y;
-20*region.y;
0*region.y;]
end
end
if true
function bcMatrix = myConstPsi4Rev(region,state)
bcMatrix = [ 0*region.y;
-20*dl*ones(size(region.y));
0*region.y ];
end
end
But the problem is because the geometry is consisting of many squares instead of one rectangle. My code becomes very slow. How can I change the geometry after meshing to the rectangle with out having any issues. Or is there any other way to do this. I can not use legacy work flow because the elements are only linear in legacy, and I need quadratic elements. My results are much smother with this structured mesh, rather than the Matlab default unstructured mesh.So I really need this structured mesh, but because of the geometry the code is super slow. Please help me.
  1 个评论
Kaveh Gharibi
Kaveh Gharibi 2017-11-29
The soloution to this problem is here:
https://www.mathworks.com/matlabcentral/answers/346410-is-there-a-way-to-import-2d-meshes-nodes-and-triangles-into-pde-model-objects

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File 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