Bug with skinny triangles

[ajt60gaibb]

h = 1e-4;
D = ultraSEMDomain.triangle([ 0 0 ; 1 0; 1 h]); 
op = ultraSEMPDO(1, 0, 0);
rhs = -1;                            
bc = 0;                              
L = ultraSEM(D, op, rhs);           
sol = L\bc;                        
plot(sol)

Error using reshape
To RESHAPE the number of
elements must not change.
Error in
ultraSEMPatch/intersect (line
103)
flags1 =
find(all(diff(reshape(i4a,na-2,numel(ia)))
== -1));
Error in ultraSEMPatch/merge
(line 20)
[i1, i2, s1, s2, flip1, flip2,
dom] = intersect(a, b);
Error in ultraSEMDomain/build
(line 107)
                        q{k} =
                        merge(pk{:});
                        %
                        Perform
                        merge.
Error in ultraSEM/build (line
129)
            S.patches =
            build(S.domain,
            S.patches);
Error in ultraSEM/solve (line
206)
            if (
            numel(S.patches) >
            1 ), build(S); end
Error in  \  (line 184)
            [varargout{1:nargout}]
            =
            solve(varargin{:}); 

[danfortunato]

In @ultraSEMPatch/intersect.m, the node coordinates are shifted by 42 before intersecting. Changing this shift to 0 seems to fix this bug: see 63f7da88b121810d6fe64989e101b0af141cd063. @nickhale, do you remember the reason for this shift, and if it is necessary?

[danfortunato]

Actually, though changing this seems to fix the skinny triangle issue, it breaks other domains. Not sure why.

[nickhale]

You're right that the 42 shift is breaking the skinny triangle, presumably because the sides are so small that when we add 42 to them in double precision the interface points are no longer unique. (Smearing error - just did this with my undergrads..)

I can't remember exactly why the 42 was necessary (although it seems it was). I'll see if I can find a different fix to the non-skinny case which doesn't involve this hack.

Alternatively, it's odd that we're still relying on points on the boundary to annihilate the interface conditions. I.e., we're still not 'grid free'. Is there some other way we can keep track of the interfaces and their indicies?

Edit: Even the single() trick is not going to be a good idea if we want really skinny elements.

Can you send an example that doesn't work if shift = 0. Everything I try seems to work..

[danfortunato]

Yeah, I think both the shift and the single() trick will be a problem for very skinny elements.

Keeping around points on the boundary made things easier for two reasons:

• Getting the indices of the coefficients correct in intersect(), which can be subtle for merging two ultraSEMParents.
• Converting the BC data at the top level into coefficients, which requires evaluating the BC data at the boundary points of the top level and converting each boundary separately into coefficients.

If you know of a more clever way to do those operations, then we could get rid of the boundary nodes.

I'll find an example that doesn't work when shift = 0...

[danfortunato]

I think the behavior I had seen before with shift = 0 may have been a clear classes issue, as I can’t reproduce any problems for shift = 0 now...

[nickhale]

That's encouraging. Do we need to worry about triangles becoming skinny/small enough that working in single precision is an issue?

If I recall correctly, the process of using intersect to find corresponding rows during the merge was horribly slow. That's why I put a lot of work in to avoiding it in the (standard) case of rectangular elements.

We haven't spoken a lot about gridding for general domains, but presumably we'd like to grid uniformly where possible with rectangles, and use triangles/kites where things get tricky? (i.e., refining corner layers if we can't do hanging nodes or for more complicated domains)

[ajt60gaibb]

I am ok with the code only being able to do triangles of height >1e-8(!)