a b/Cobiveco/@cobiveco/prepareMesh2.m
1
function prepareMesh2(o)
2
3
if ~o.available.transmural
4
    o.computeTransmural;
5
end
6
o.printStatus('Preparing mesh 2...');
7
t = toc;
8
9
isovalue = 0.5;
10
numTries = 5;
11
for i = 1:numTries
12
    [o.m2.vol,mmgStatus,mmgOutput2] = mmg(o.m1.vol, o.m1.ridgeLaplace, sprintf('-ls %1.5e -nr -hausd %1.5e -hmin %1.5e -hmax %1.5e', isovalue, o.cfg.mmgSizingParam(:)'*o.m0.meanEdgLen));
13
14
    if o.cfg.exportLevel > 1 || mmgStatus ~= 0
15
        if i == 1
16
            fid = fopen(sprintf('%smmgOutput2.txt', o.cfg.outPrefix), 'w');
17
        else
18
            fid = fopen(sprintf('%smmgOutput2.txt', o.cfg.outPrefix), 'a');
19
        end
20
        fprintf(fid, '%s', mmgOutput2);
21
        fclose(fid);
22
    end
23
24
    if mmgStatus == 0
25
        break;
26
    elseif i < numTries
27
        warning('Mmg remeshing with an isovalue of %.3f failed (system command status %i). Trying again with a slightly larger isovalue.', isovalue, mmgStatus);
28
        isovalue = isovalue + 0.1*prctile(max(o.m1.ridgeLaplace(o.m1.vol.cells),[],2)-min(o.m1.ridgeLaplace(o.m1.vol.cells),[],2),95);
29
        o.cfg.mmgSizingParam(1) = 0.8*o.cfg.mmgSizingParam(1);
30
    else
31
        error('Mmg remeshing failed ultimately after trying %i different isovalues (system command status %i).', numTries, mmgStatus);
32
    end
33
end
34
35
o.m2.sur = vtkDataSetSurfaceFilter(vtkDeleteDataArrays(o.m2.vol));
36
o.m2.sur = vtkArrayMapperNearestNeighbor(o.m0.sur, o.m2.sur);
37
o.m2.surToVol = vtkMapPointIds(o.m2.vol, o.m2.sur);
38
39
P2 = double(o.m2.vol.points);
40
C2 = double(o.m2.vol.cells);
41
o.m2.L = cotmatrix(P2, C2);
42
o.m2.G = grad(P2, C2);
43
o.m2.M = baryInterpMat(P2, C2, o.m0.vol.points);
44
45
if o.cfg.exportLevel > 1
46
    o.m2.debug = o.m2.vol;
47
    o.m2.debug.pointData.surClass = repmat(uint8(0),size(o.m2.debug.points,1),1);
48
    o.m2.debug.pointData.surClass(o.m2.surToVol) = o.m2.sur.pointData.class;
49
end
50
51
o.printStatus(sprintf('%.1f seconds\n', toc-t), true);
52
o.available.mesh2 = true;
53
54
end