Merge branch 'master' into documentation-paper

lib/assembly/@BilinearForm/BilinearForm.m

@@ -4,11 +4,11 @@
 
 classdef BilinearForm < handle    
     %% properties
-    properties (GetAccess='public', SetAccess='protected')
+    properties (GetAccess=public, SetAccess=protected)
         fes
     end
     
-    properties (Access='public')
+    properties (Access=public)
         a {mustBeEvaluableOrEmpty} = []     % Diffusion matrix (or scalar)
         b {mustBeEvaluableOrEmpty} = []     % Convection vector field
         c {mustBeEvaluableOrEmpty} = []     % Reaction coefficient
@@ -33,11 +33,11 @@ classdef BilinearForm < handle
         end
     end
     
-    methods (Access='public')
+    methods (Access=public)
         mat = assemble(obj);
     end
     
-    methods (Static, Access='protected')
+    methods (Static, Access=protected)
         val = diffusionPart(a, Dphi);
         val = convectionPart(b, Dphi, phi);
         val = reactionPart(c, phi);

lib/assembly/@LinearForm/LinearForm.m

@@ -4,11 +4,11 @@
 
 classdef LinearForm < handle
     %% properties
-    properties (GetAccess='public', SetAccess='protected')
+    properties (GetAccess=public, SetAccess=protected)
         fes
     end
     
-    properties (Access='public')
+    properties (Access=public)
         f {mustBeEvaluableOrEmpty} = []         % Volume load
         fvec {mustBeEvaluableOrEmpty} = []      % Generalized volume load
         robin {mustBeEvaluableOrEmpty} = []     % Robin boundary load
@@ -33,11 +33,11 @@ classdef LinearForm < handle
         end
     end
     
-    methods (Access='public')
+    methods (Access=public)
         vec = assemble(obj);
     end
     
-    methods (Static, Access='protected')
+    methods (Static, Access=protected)
         val = scalarPart(f, phi);
         val = vectorPart(fvec, Dphi);
         val = neumannPart(neumann, phi);

lib/functions/@Evaluable/Evaluable.m

@@ -8,16 +8,16 @@
 
 classdef Evaluable < handle
     %% properties
-    properties (Abstract, GetAccess='public', SetAccess='protected')
+    properties (Abstract, GetAccess=public, SetAccess=protected)
         mesh
     end
     
     %% methods
-    methods (Abstract, Access='public')
+    methods (Abstract, Access=public)
         eval(obj, bary, idx)
     end
     
-    methods (Access='public')
+    methods (Access=public)
         plot(obj, plotOptions)
         evalEdge(obj, bary, idx)
     end

lib/functions/@Evaluable/plot.m

@@ -59,18 +59,17 @@ else
         warning('Custom plot options discarded due to subdivisions.')
     end
     for s = smin:smax
-        idx = find(subdivision == s);
-        if ~isempty(idx)
+        currentElems = find(subdivision == s);
+        if ~isempty(currentElems)
             % subtriangles on one element
             bary = Barycentric2D(partition3(s)/s);
             [elements, outerEdges] = getSubtriangles(s);
-            nSubtriangles = size(elements, 1) * numel(idx);
             
             % compute plot data
-            coordinates = mesh.elementwiseCoordinates(bary, idx);
+            coordinates = mesh.elementwiseCoordinates(bary, currentElems);
             Xel = splitFirstDimension(coordinates, elements);
             Xlin = splitFirstDimension(coordinates, outerEdges);
-            data = guaranteeSize(eval(obj, bary, idx), mesh.nElements, bary.nNodes);
+            data = guaranteeSize(eval(obj, bary, currentElems), numel(currentElems), bary.nNodes);
             Zel = splitFirstDimension(data, elements);
             Zlin = splitFirstDimension(data, outerEdges);
             

lib/functions/CompositeFunction.m

@@ -13,14 +13,14 @@
 
 classdef CompositeFunction < Evaluable
     %% properties
-    properties (SetAccess='protected', GetAccess='public')
+    properties (SetAccess=protected, GetAccess=public)
         mesh
         functionHandle
         functionArguments
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = CompositeFunction(functionHandle, functionArguments)
             arguments
                 functionHandle (1,1) function_handle
@@ -44,7 +44,7 @@ classdef CompositeFunction < Evaluable
         end
     end
     
-    methods (Access='private')
+    methods (Access=private)
         function val = compositeEval(obj, evalMethod, bary, idx)
             arguments
                 obj

lib/functions/Constant.m

@@ -6,13 +6,13 @@
 
 classdef Constant < Evaluable
     %% properties
-    properties (GetAccess='public', SetAccess='protected')
+    properties (GetAccess=public, SetAccess=protected)
         mesh
         constantValue (:,1) double
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = Constant(mesh, value)
             obj.mesh = mesh;
             obj.constantValue = value;

lib/functions/FeFunction.m

@@ -19,19 +19,19 @@
 
 classdef FeFunction < Evaluable
     %% properties
-    properties (SetAccess='protected', GetAccess='public')
+    properties (SetAccess=protected, GetAccess=public)
         mesh
         fes
         data (1,:) double
         dataIsValid (1,1) logical
     end
     
-    properties (Access = 'private')
+    properties (Access = private)
         listenerHandle
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = FeFunction(fes)
             arguments
                 fes (1,1) FeSpace
@@ -43,7 +43,7 @@ classdef FeFunction < Evaluable
             obj.dataIsValid = true;
             
             % register FeFunction as an observer of mesh
-            obj.listenerHandle = fes.mesh.listener('HasChanged', @obj.handleChangedMesh);
+            obj.listenerHandle = fes.mesh.listener('JustRefined', @obj.handleChangedMesh);
         end
         
         function setData(obj, data)
@@ -91,7 +91,7 @@ classdef FeFunction < Evaluable
         end
     end
     
-    methods (Access = 'private')
+    methods (Access = private)
         function handleChangedMesh(obj, ~, ~)
             obj.dataIsValid = false;
         end

lib/functions/Gradient.m

@@ -8,13 +8,13 @@
 
 classdef Gradient < Evaluable
     %% properties
-    properties (SetAccess='protected', GetAccess='public')
+    properties (SetAccess=protected, GetAccess=public)
         u
         mesh
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = Gradient(u)
             arguments
                 u (1,1) FeFunction

lib/functions/Hessian.m

@@ -8,13 +8,13 @@
 
 classdef Hessian < Evaluable
     %% properties
-    properties (SetAccess='protected', GetAccess='public')
+    properties (SetAccess=protected, GetAccess=public)
         u
         mesh
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = Hessian(u)
             arguments
                 u (1,1) FeFunction

lib/functions/MeshFunction.m

@@ -8,13 +8,13 @@
 
 classdef MeshFunction < Evaluable
     %% properties
-    properties (GetAccess='public', SetAccess='protected')
+    properties (GetAccess=public, SetAccess=protected)
         mesh
         functionHandle
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = MeshFunction(mesh, functionHandle)
             arguments
                 mesh (1,1) Mesh

lib/functions/NormalVector.m

@@ -8,12 +8,12 @@
 
 classdef NormalVector < Evaluable
     %% properties
-    properties (GetAccess='public', SetAccess='protected')
+    properties (GetAccess=public, SetAccess=protected)
         mesh
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = NormalVector(mesh)
             obj.mesh = mesh;
         end

lib/functions/TestFunction.m

@@ -9,13 +9,13 @@
 
 classdef TestFunction < Evaluable
     %% properties
-    properties (SetAccess='protected', GetAccess='public')
+    properties (SetAccess=protected, GetAccess=public)
         fes
         mesh
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = TestFunction(fes)
             arguments
                 fes (1,1) FeSpace

lib/functions/TestFunctionGradient.m

@@ -10,13 +10,13 @@
 
 classdef TestFunctionGradient < Evaluable
     %% properties
-    properties (SetAccess='protected', GetAccess='public')
+    properties (SetAccess=protected, GetAccess=public)
         fes
         mesh
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = TestFunctionGradient(fes)
             arguments
                 fes (1,1) FeSpace

lib/mesh/@Mesh/Mesh.m

@@ -6,7 +6,7 @@
 
 classdef Mesh < handle
     %% properties
-    properties (GetAccess = 'public', SetAccess = 'private')
+    properties (GetAccess = public, SetAccess = private)
         coordinates (2,:) double
         elements (3,:) double
         edges (2,:) double
@@ -15,11 +15,6 @@ classdef Mesh < handle
         boundaries (:,1) cell
     end
     
-    properties (GetAccess = 'public', SetAccess = 'private', Hidden = true)
-        markedEdges (:,1) logical
-        bisecGroups (:,1) cell
-    end
-    
     %% dependent properties (computed from data)
     properties (Dependent)
         nCoordinates
@@ -29,18 +24,19 @@ classdef Mesh < handle
     end
     
     %% private properties for caching purposes
-    properties (Access = 'private')
+    properties (Access = private)
         trafo
     end
     
     %% events
     events
         IsAboutToRefine
-        HasChanged
+        JustRefined
+        RefineCompleted
     end
     
     %% public methods
-    methods (Access = 'public')
+    methods (Access = public)
         function obj = Mesh(coordinates, elements, boundaries)
             % Construct Mesh object from given coordinate, element, and boundary
             % arrays.
@@ -70,12 +66,13 @@ classdef Mesh < handle
         saveTikzConforming(obj, varargin)
         refineLocally(obj, marked, method)
         refineUniform(obj, n, method)
+        changeRefinementEdge(obj, newRefinementEdge)
         edges = getCombinedBndEdges(obj, idx)
     end
     
     %% protected methods
-    methods (Access='protected')
-        updateData(obj)
+    methods (Access=protected)
+        updateData(obj, bisecData)
     end
     
     %% get methods for dependent data

lib/mesh/@Mesh/changeRefinementEdge.m

+% changeRefinementEdge Rotate internal data structures such that given
+%   edge is new refinement edge per element.
+%
+% mesh.changeRefinementEdge(newRefinementEdge) changes the elements,
+%   element2edges, and flipEdges data structures of the mesh such that the
+%   new edge is the edge given in newRefinementEdge per element. New edges
+%   are given as index in {1,2,3}, according to the old element2edges data
+%   structure.
+
+function changeRefinementEdge(obj, newRefinementEdge)
+
+arguments
+    obj
+    newRefinementEdge (1,:) double {mustBeInRange(newRefinementEdge, 1, 3)}
+end
+
+for k = 2:3
+    idx = find(newRefinementEdge == k);
+    obj.element2edges(:,idx) = circshift(obj.element2edges(:,idx), 4-k, 1);
+    obj.elements(:,idx) = circshift(obj.elements(:,idx), 4-k, 1);
+    obj.flipEdges(:,idx) = circshift(obj.flipEdges(:,idx), 4-k, 1);
+end
+
+end
\ No newline at end of file

lib/mesh/@Mesh/refineLocally.m

@@ -22,14 +22,14 @@ end
 
 %% prepare data for refinement and make it available to listeners
 refinementAlgorithm = feval(method, obj);
-[obj.markedEdges, obj.bisecGroups] = refinementAlgorithm.prepareRefinementData(marked);
-obj.notify('IsAboutToRefine')
+bisecData = refinementAlgorithm.prepareRefinementData(marked);
+obj.notify('IsAboutToRefine', bisecData)
 
 %%  do actual refinement and clear data of old mesh
-obj.updateData();
-obj.markedEdges = [];
-obj.bisecGroups = [];
+obj.updateData(bisecData);
 obj.trafo = [];
-obj.notify('HasChanged');
+obj.notify('JustRefined');
+
+obj.notify('RefineCompleted')
 
 end
\ No newline at end of file

lib/mesh/@Mesh/updateData.m

@@ -5,80 +5,86 @@
 %
 %   See also Mesh
 
-function updateData(obj)
-%% create indices of ...
-% new nodes
-edgeNodes = obj.nCoordinates + (1:nnz(obj.markedEdges))';
-innerNodes = cumsum([edgeNodes(end); cellfun(@(x) x.nMembers*x.nInnerNodes, obj.bisecGroups)]);
-
-% new edges
-innerEdges = cumsum([obj.nEdges + nnz(obj.markedEdges); cellfun(@(x) x.nMembers*x.nInnerEdges, obj.bisecGroups)]);
-nChildEdges = 1 + obj.markedEdges;
-oldEdge2parent = cumsum([1; nChildEdges]);
-
-% new elements
-nChildElems = ones(obj.nElements, 1);
-for k = 1:length(obj.bisecGroups)
-    nChildElems(obj.bisecGroups{k}.elementIdx) = obj.bisecGroups{k}.nDescendants;
-end
-oldElement2parent = cumsum([1; nChildElems]);
+function updateData(obj, bisecData)
+%% create old-to-new connectivity
+oldEdge2parent = cumsum([1; 1+bisecData.bisectedEdges]);
+oldElement2parent = cumsum([1; getNChildrenPerElement(bisecData)]);
 
 %% create new coordinates
-newCoordinates = [obj.coordinates, ...
-    edgewiseCoordinates(obj, [1;1]/2, obj.markedEdges), ...
-    zeros(2, innerNodes(end)-innerNodes(1))];
-
-for k = 1:numel(obj.bisecGroups)
-    if obj.bisecGroups{k}.nInnerNodes ~= 0
-        newCoordinates(:,innerNodes(k)+1:innerNodes(k+1)) = ...
-            elementwiseCoordinates(obj, obj.bisecGroups{k}.innerNodes, obj.bisecGroups{k}.elementIdx);
-    end
+% new nodes = old nodes | nodes on edges | interior nodes
+idxN = cumsum([1; obj.nCoordinates; nnz(bisecData.bisectedEdges); ...
+    bisecData.nRefinedElements.*bisecData.nInnerNodes]);
+
+newCoordinates = zeros(2, idxN(end)-1);
+newCoordinates(:,idxN(1):idxN(2)-1) = obj.coordinates;
+newCoordinates(:,idxN(2):idxN(3)-1) = edgewiseCoordinates(obj, [1;1]/2, bisecData.bisectedEdges);
+
+for k = find(bisecData.nInnerNodes ~= 0)'
+    newCoordinates(:,idxN(k+2):idxN(k+3)-1) = ...
+        elementwiseCoordinates(obj, bisecData.bisection{k}.innerNodes, getRefinedElementIdx(bisecData, k));
 end
 
 %% refine existing edges
-newEdges = [moveParents(obj.edges, oldEdge2parent), ...
-    zeros(2, innerEdges(end) - innerEdges(1))];
+% new edges = old edges + children | interior edges
+idxS = cumsum([1; obj.nEdges + nnz(bisecData.bisectedEdges); ...
+    bisecData.nRefinedElements.*bisecData.nInnerEdges]);
 
-children = getChildIndices(oldEdge2parent(obj.markedEdges), 2);
-newEdges(:,children) = [obj.edges(1,obj.markedEdges), edgeNodes'; ...
-    edgeNodes', obj.edges(2,obj.markedEdges)];
+newEdges = zeros(2, idxS(end)-1);
+newEdges(:,oldEdge2parent(1:(end-1))) = obj.edges;
+
+edgeNodes = (idxN(2):idxN(3)-1)';
+children = getChildIndices(oldEdge2parent(bisecData.bisectedEdges), 2);
+newEdges(:,children) = [obj.edges(1,bisecData.bisectedEdges), edgeNodes'; ...
+    edgeNodes', obj.edges(2,bisecData.bisectedEdges)];
 
 %% refine boundaries
 newBoundaries = cell(obj.nBoundaries, 1);
 for k = 1:obj.nBoundaries
     newBoundaries{k} = oldEdge2parent(obj.boundaries{k});
-    bndChildren = newBoundaries{k}(obj.markedEdges(obj.boundaries{k})) + 1;
+    bndChildren = newBoundaries{k}(bisecData.bisectedEdges(obj.boundaries{k})) + 1;
     newBoundaries{k} = [newBoundaries{k}; bndChildren];
 end
 
-%% allocate new element-wise arrays
-oldEdgeNewIdx = oldEdge2parent(obj.element2edges);
-newElements = moveParents(obj.elements, oldElement2parent);
-newE2E = moveParents(oldEdgeNewIdx, oldElement2parent);
-newFlip = moveParents(obj.flipEdges, oldElement2parent);
+%% refine element data
+newElements = zeros(3, oldElement2parent(end)-1);
+newE2E = zeros(3, oldElement2parent(end)-1);
+newFlip = false(3, oldElement2parent(end)-1);
 
-%% generate bisection rule dependent data
-[left, right] = getLeftAndRightEdgePerElement(obj, oldEdgeNewIdx);
+oldEdgeNewIdx = oldEdge2parent(obj.element2edges);
+elems = bisecData.getRefinedElementIdx(0);
+newElements(:,oldElement2parent(elems)) = obj.elements(:,elems);
+newE2E(:,oldElement2parent(elems)) = oldEdgeNewIdx(:,elems);
+newFlip(:,oldElement2parent(elems)) = obj.flipEdges(:,elems);
+
+%% get left and right edge per element
+% here, we use the fact that left and right edge index only differ by 1
+bisectedEdges = bisecData.bisectedEdges(obj.element2edges);
+left = oldEdgeNewIdx.*bisectedEdges;
+right = left + bisectedEdges;
+left = left + obj.flipEdges;
+right = right - obj.flipEdges;
+
+%% refine elements + add inner edges
 edgeMidPts = zeros(obj.nEdges,1);
-edgeMidPts(obj.markedEdges) = edgeNodes;
+edgeMidPts(bisecData.bisectedEdges) = edgeNodes;
 edgeMidPts = edgeMidPts(obj.element2edges);
 
-for k = 1:length(obj.bisecGroups)
-    idx = (innerEdges(k)+1):innerEdges(k+1);
-    intEdges = reshape(idx, [], obj.bisecGroups{k}.nInnerEdges)';
-    intNodes = reshape((innerNodes(k)+1:innerNodes(k+1))', obj.bisecGroups{k}.nInnerNodes, []);
+for k = find(bisecData.nRefinedElements)'
+    idx = idxS(k+1):idxS(k+2)-1;
+    elems = getRefinedElementIdx(bisecData, k);
+    intEdges = reshape(idx, [], bisecData.nInnerEdges(k))';
+    intNodes = reshape((idxN(k+2):idxN(k+3)-1)', bisecData.nInnerNodes(k), []);
     
-    parentEdges = oldEdgeNewIdx(:, obj.bisecGroups{k}.elementIdx);
-    parentElements = oldElement2parent(obj.bisecGroups{k}.elementIdx);
-    children = getChildIndices(parentElements, obj.bisecGroups{k}.nDescendants);
+    parentEdges = oldEdgeNewIdx(:, elems);
+    parentElements = oldElement2parent(elems);
+    children = getChildIndices(parentElements, bisecData.nDescendants(k));
     
-    arg = restrictTo(obj.bisecGroups{k}.elementIdx, ...
-        obj.elements, edgeMidPts, left, right, obj.flipEdges);
+    arg = restrictTo(elems, obj.elements, edgeMidPts, left, right, obj.flipEdges);
     
-    newEdges(:,idx) = obj.bisecGroups{k}.createNewEdges(arg{1}, arg{2}, intNodes);
-    newElements(:,children) = obj.bisecGroups{k}.refineElement(arg{1}, arg{2}, intNodes);
-    newE2E(:,children) = obj.bisecGroups{k}.refineEdgeConnectivity(parentEdges, arg{3}, arg{4}, intEdges);
-    newFlip(:,children) = obj.bisecGroups{k}.refineEdgeOrientation(arg{5});
+    newEdges(:,idx) = bisecData.bisection{k}.createNewEdges(arg{1}, arg{2}, intNodes);
+    newElements(:,children) = bisecData.bisection{k}.refineElement(arg{1}, arg{2}, intNodes);
+    newE2E(:,children) = bisecData.bisection{k}.refineEdgeConnectivity(parentEdges, arg{3}, arg{4}, intEdges);
+    newFlip(:,children) = bisecData.bisection{k}.refineEdgeOrientation(arg{5});
 end
 
 %% assign new arrays
@@ -96,20 +102,6 @@ function children = getChildIndices(parents, nDescendants)
     children = reshape((0:(nDescendants-1)) + parents, [], 1);
 end
 
-function newArray = moveParents(oldArray, parentPositions)
-    newArray = zeros(size(oldArray, Dim.Vector), parentPositions(end)-1, class(oldArray));
-    newArray(:,parentPositions(1:(end-1))) = oldArray;
-end
-
 function restrictedArrays = restrictTo(idx, varargin)
     restrictedArrays = cellfun(@(x) x(:,idx), varargin, 'UniformOutput', false);
-end
-
+end
\ No newline at end of file
-function [left, right] = getLeftAndRightEdgePerElement(mesh, oldEdgeNewIdx)
-    % here, we use the fact that left and right edge index only differ by 1
-    bisectedEdges = mesh.markedEdges(mesh.element2edges);
-    left = oldEdgeNewIdx.*bisectedEdges;
-    right = left + bisectedEdges;
-    left = left + mesh.flipEdges;
-    right = right - mesh.flipEdges;
-end

lib/mesh/refinements/NVB.m

@@ -7,24 +7,24 @@
 
 classdef NVB
     %% properties
-    properties (GetAccess='public', SetAccess='protected')
+    properties (GetAccess=public, SetAccess=protected)
         mesh
     end
     
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = NVB(mesh)
             obj.mesh = mesh;
         end
         
-        function [markedEdges, bisecGroups] = prepareRefinementData(obj, markedElements)
+        function bisecData = prepareRefinementData(obj, markedElements)
             markedEdges = obj.markedElementsToEdges(markedElements);
             markedEdges = obj.meshClosure(markedEdges);
-            bisecGroups = obj.groupElements(markedEdges, markedElements);
+            bisecData = obj.groupElements(markedEdges, markedElements);
         end
     end
     
-    methods (Access='protected')
+    methods (Access=protected)
         function markedEdges = markedElementsToEdges(obj, markedElements)
             markedEdges = false(obj.mesh.nEdges,1);
             markedEdges(obj.mesh.element2edges(:,markedElements)) = true;
@@ -39,13 +39,14 @@ classdef NVB
             end
         end
         
-        function bisecGroups = groupElements(obj, markedEdges, ~)
+        function bisecData = groupElements(obj, markedEdges, ~)
             % group elements according to their refinement pattern
             edge = reshape(markedEdges(obj.mesh.element2edges), size(obj.mesh.element2edges));
-            bisecGroups{1} = Bisec1  ( find( edge(1,:) & ~edge(2,:) & ~edge(3,:) ) );
-            bisecGroups{2} = Bisec12 ( find( edge(1,:) &  edge(2,:) & ~edge(3,:) ) );
-            bisecGroups{3} = Bisec13 ( find( edge(1,:) & ~edge(2,:) &  edge(3,:) ) );
-            bisecGroups{4} = Bisec123( find( edge(1,:) &  edge(2,:) &  edge(3,:) ) );
+            bisecData = BisectionEventData(obj.mesh.nElements, markedEdges, ...
+                Bisec1,   find(edge(1,:) & ~edge(2,:) & ~edge(3,:)), ...
+            	Bisec12,  find(edge(1,:) &  edge(2,:) & ~edge(3,:)), ...
+            	Bisec13,  find(edge(1,:) & ~edge(2,:) &  edge(3,:)), ...
+            	Bisec123, find(edge(1,:) &  edge(2,:) &  edge(3,:)));
         end
     end
 end

lib/mesh/refinements/NVB1.m

@@ -3,13 +3,13 @@
 
 classdef NVB1 < NVB
     %% methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = NVB1(mesh)
             obj = obj@NVB(mesh);
         end
     end
     
-    methods (Access='protected')
+    methods (Access=protected)
         function markedEdges = markedElementsToEdges(obj, markedElements)
             markedEdges = false(1,obj.mesh.nEdges);
             markedEdges(obj.mesh.element2edges(1,markedElements)) = true;

lib/mesh/refinements/NVB3.m

@@ -3,9 +3,9 @@
 
 classdef NVB3 < NVB   
     %% public methods
-    methods (Access='public')
+    methods (Access=public)
         function obj = NVB3(mesh)
             obj = obj@NVB(mesh);
         end
     end
-end
+end
\ No newline at end of file