Clear Filters
Clear Filters

Pde toolbox: internal decsg error

2 views (last 30 days)
Brian
Brian on 4 Jan 2015
Edited: Brian on 4 Jan 2015
Hello,
I'm using the decsg function to create a mesh with random circular inclusions. I have to do Monte Carlo simulations in order to estimate some mean values and standard deviations. It works fine but sometimes I get an "internal decsg error, (line 1435)". It occurs randomly and it is quite annoying since I have to generate at least 10000 realizations of the mesh.
Here is my code for generating the mesh with decsg:
function [p, t, ConnectElemPhase, BoundaryNodes] = ...
ConstructMesh(rayon, L, nbi, hmax)
%nbi: number of inclusions
%rayon: radius of the inclusions
%L: length of the representative volume element
a = -L+rayon+0.01; b = L-rayon-0.01;
Cadre = [3, 4, -L, L, L, -L, -L, -L, L, L]';
centre = zeros(nbi, 2);
inclusions = zeros(10, nbi);
centre(1, :) = a + (b-a)*rand(2, 1); %center of the first inclusion
inclusions(:, 1) = [1, centre(1, :), rayon, 0, 0, 0, 0, 0, 0]';
SF = 'R';
C(1, 1) = {strcat(['C', num2str(1)])};
Ci = strcat(['+C', num2str(1)]);
SF = strcat([SF, num2str(Ci)]);
for i=2:nbi
centre(i, :) = a + (b-a)*rand(2, 1); %center of the i-th inclusion
dr = sqrt((centre(1:i-1, 1)-centre(i, 1)).^2 + ...
(centre(1:i-1, 2)-centre(i, 2)).^2);
%dr: distance between the i-th inclusion and the other inclusions 1, ..., i-1.
while any(dr<=2.1*rayon)
% avoid intersection between the inclusions
% if intersection, generate a new center
centre(i, :) = a + (b-a)*rand(2, 1);
dr = sqrt((centre(1:i-1, 1)-centre(i, 1)).^2 + ...
(centre(1:i-1, 2)-centre(i, 2)).^2);
end
inclusions(:, i) = [1, centre(i, :), rayon, 0, 0, 0, 0, 0, 0]';
C(1, i) = {strcat(['C', num2str(i)])};
Ci = strcat(['+C', num2str(i)]);
SF = strcat([SF, num2str(Ci)]);
end
NS = [char('R', char(C))]'
GD = [Cadre, inclusions]; g = decsg(GD, SF, NS);
[p, e, t] = initmesh(g, 'Hmax', hmax, 'Hgrad', 1.2);
%...rest of the code not useful for my problem
Thank you for your help!
Brian.

Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!