Make access to properties protected (accessible to children).

+mgen/MgGeoArea4Line.m

 classdef MgGeoArea4Line < mgen.MgGeoArea
     % A 4-line defined area of the geometry
 
-    properties
+    properties (Access = protected)
         geoLines  = mgen.MgGeoLine.empty(0,4);
 
         gridNodes = mgen.MgGridNode.empty(0,0);
@@ -67,45 +67,51 @@ classdef MgGeoArea4Line < mgen.MgGeoArea
             self.geoLines(4).switchStartGridNode();
 
             eps = 1.e-6;
-tic
+            
             % create grid nodes
-            for ii = 2:dim0-1
+            for ii = 2:dim0-1 % straight edges
+
+                % determine local coordinates of "same" points on both straight edges
+                A = self.geoLines(1).getGridNode(   1).coord;
+                B = self.geoLines(1).getGridNode(dim0).coord;
+                C = self.geoLines(1).getGridNode(  ii).coord;
+                
+                z1 = norm((C-A)/(B-A));
 
-                N = zeros(1,2);
-                a = 0; b = 0; c = 0; d = 0;
+                A = self.geoLines(3).getGridNode(   1).coord;
+                B = self.geoLines(3).getGridNode(dim0).coord;
+                C = self.geoLines(3).getGridNode(  ii).coord;
+                
+                z3 = norm((C-A)/(B-A));
                 
                 pos = (ii-1)*dim1 + 2;
                 
-                A  = self.geoLines(1).getGridNode(ii).coord;
-                AB = self.geoLines(3).getGridNode(ii).coord - A;
+                for jj = 2:dim1-1 % curved edges
 
-                if abs(AB(1)) <= eps
-                    N(1) = A(1);
-                else
-                    a = AB(2) / AB(1);
-                    b = A(2) - a * A(1);
-                end
+                    % determine local coordinates of "same" points on both straight edges
+                    A = self.geoLines(2).getGridNode(   1).coord;
+                    B = self.geoLines(2).getGridNode(dim1).coord;
+                    C = self.geoLines(2).getGridNode(  jj).coord;
 
-                for jj = 2:dim1-1
+                    z2 = norm((C-A)/(B-A));
 
-                    C  = self.geoLines(2).getGridNode(jj).coord;
-                    CD = self.geoLines(4).getGridNode(jj).coord - C;
+                    A = self.geoLines(4).getGridNode(   1).coord;
+                    B = self.geoLines(4).getGridNode(dim1).coord;
+                    C = self.geoLines(4).getGridNode(  jj).coord;
 
-                    if abs(CD(1)) <= eps
-                        N(1) = C(1);
-                    else
-                        c = CD(2) / CD(1);
-                        d = C(2) - c * C(1);
-                    end
+                    z4 = norm((C-A)/(B-A));
                     
-                    if abs(AB(1)) <= eps && abs(CD(1)) > eps
-                        N(2) = c * N(1) + d;
-                    elseif abs(CD(1)) <= eps && abs(AB(1)) > eps
-                        N(2) = a * N(1) + b;
-                    else
-                        N(1) = (d-b) / (a-c);
-                        N(2) = a * N(1) + b;
-                    end
+                    A = self.geoLines(2).getGridNode(jj).coord;
+                    B = self.geoLines(4).getGridNode(jj).coord;
+                    C = self.geoLines(1).getGridNode(ii).coord;
+                    D = self.geoLines(3).getGridNode(ii).coord;
+                    
+                    zi = (ii-1)/(dim0-1); % local z along curve
+                    zj = (jj-1)/(dim1-1); % local z along curve
+                    zzi = (1-zi) * z2 + zi * z4; % interpolate outer pos
+                    zzj = (1-zj) * z1 + zj * z3; 
+                    
+                    N = ( (1-zzj)*A + (zzj)*B + (1-zzi)*C + (zzi)*D )/2;
 
                     % position node
                     self.gridNodes(pos) = mgen.MgGridNode(run_id, N);
@@ -116,7 +122,7 @@ tic
                     pos = pos + 1;
                 end
             end
-toc
+            
             % sort in existing nodes
             for ii = 1:dim1
                 self.gridNodes(ii) = self.geoLines(2).getGridNode(ii);