diff --git a/examples/goalOrientedAFEM.m b/examples/goalOrientedAFEM.m
index c270eba74963aad09d60c3e803885ae2a59e6af8..e966ffa0bc7ff92cf2f609649d2436c66ea0e58c 100644
--- a/examples/goalOrientedAFEM.m
+++ b/examples/goalOrientedAFEM.m
@@ -39,8 +39,7 @@ for p = 1:pmax
lfG.qrfvec = QuadratureRule.ofOrder(max(p-1, 1));
%% set up lifting operators for rhs FEM-data
- Pv = LoL2Prolongation(v);
- Pw = LoL2Prolongation(w);
+ P = LoFeProlongation(ncFes);
%% adaptive loop
i = 1;
@@ -70,8 +69,8 @@ for p = 1:pmax
marked = markGoafemMS(eta2, zeta2, theta);
mesh.refineLocally(marked, 'NVB');
i = i+1;
- v.setData(Pv.interpolatedData);
- w.setData(Pw.interpolatedData);
+ v.setData(prolongate(P, v));
+ w.setData(prolongate(P, w));
end
end
diff --git a/examples/iterativeLinearization.m b/examples/iterativeLinearization.m
index 7e199187ec144f0b0ec32084bb59f5c9ea6072bb..9c9ea5efaf0902140faec7b93cc08eefcb30ef3e 100644
--- a/examples/iterativeLinearization.m
+++ b/examples/iterativeLinearization.m
@@ -24,7 +24,10 @@ for k = 1:length(linearizations)
g = Constant(mesh, 1);
[blf, lf] = setProblemData(mesh, fes, Du, g, linearizations(k));
eDensity = CompositeFunction(@(p) muIntegral(vectorProduct(p, p)), Du);
- u.setData(0);
+ u.setData(0);
+
+ %% nested iteration
+ P = LoFeProlongation(fes);
%% adaptive algorithm
i = 1;
@@ -67,13 +70,10 @@ for k = 1:length(linearizations)
break
end
- %% nested iteration
- uold = LoH1Prolongation(u);
-
%% refine mesh
marked = markDoerflerSorting(eta2, theta);
mesh.refineLocally(marked, 'NVB');
- u.setData(uold.interpolatedData);
+ u.setData(prolongate(P, u));
i = i+1;
end
end
diff --git a/lib/assembly/@BilinearForm/BilinearForm.m b/lib/assembly/@BilinearForm/BilinearForm.m
index 866ad7269fd88aa6d68171b2ea6c34f87fe0dc7d..fbd4f0e807b96750ef9c6f8bf92fb44a780a6f50 100644
--- a/lib/assembly/@BilinearForm/BilinearForm.m
+++ b/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);
diff --git a/lib/assembly/@LinearForm/LinearForm.m b/lib/assembly/@LinearForm/LinearForm.m
index c79c0efd22c6a3066a5d6045d7cc24f5f4c909d0..7fc9ddc830b59b8297d949a37b77d16d9a0d2b9b 100644
--- a/lib/assembly/@LinearForm/LinearForm.m
+++ b/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);
diff --git a/lib/functions/@Evaluable/Evaluable.m b/lib/functions/@Evaluable/Evaluable.m
index f42571dad89a4d5fec6899694c6e1d80504a2ad7..683d3516f77e94165b7c7c7e454f41d3b0d1e613 100644
--- a/lib/functions/@Evaluable/Evaluable.m
+++ b/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
diff --git a/lib/functions/CompositeFunction.m b/lib/functions/CompositeFunction.m
index a62ed3e0eb396a640eef55a6eefbcc3a137c49ef..5b5ef2f5160c390acb6df31e22cf98803d44e2a7 100644
--- a/lib/functions/CompositeFunction.m
+++ b/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
diff --git a/lib/functions/Constant.m b/lib/functions/Constant.m
index 324a00750020ca5f8fcd865567d834bc1798b74a..28c902afe17aa69fc4c3761f0cd494e93147864c 100644
--- a/lib/functions/Constant.m
+++ b/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;
diff --git a/lib/functions/FeFunction.m b/lib/functions/FeFunction.m
index d47d0c91a7aa2e4a6ae4cc25e3662b7070330f11..510a0239132f4f1c69ecf5776fb2a47999cf320b 100644
--- a/lib/functions/FeFunction.m
+++ b/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
diff --git a/lib/functions/Gradient.m b/lib/functions/Gradient.m
index 8eae3d8357beee37b70f6c29f8ff10d5e734d291..46b34ae96ff394dabc3e4004d2297fa99a8fcd33 100644
--- a/lib/functions/Gradient.m
+++ b/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
diff --git a/lib/functions/Hessian.m b/lib/functions/Hessian.m
index 74d69d22dfcf1cf761894189e48f3447e2d4f1de..35a75fc2ee6cd2ba20aa6107c3fb1b5e806ac1f9 100644
--- a/lib/functions/Hessian.m
+++ b/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
diff --git a/lib/functions/MeshFunction.m b/lib/functions/MeshFunction.m
index 6f52fb3a2ec49fc78732a49a9c81f3a3c97a311d..df0d0028cba8914c677724479d8c419eefe0c874 100644
--- a/lib/functions/MeshFunction.m
+++ b/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
diff --git a/lib/functions/NormalVector.m b/lib/functions/NormalVector.m
index c19d24a1ded1d5454840da9578413a18a57d4cdc..2c728cc7bdff59a440139ef69ab5abe76865d278 100644
--- a/lib/functions/NormalVector.m
+++ b/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
diff --git a/lib/functions/TestFunction.m b/lib/functions/TestFunction.m
index 7e2cbeeea44f4a6f7486a475fc37b1e2b1b6bbae..15a636c0fd1ff85f22cac60fb209530956b14540 100644
--- a/lib/functions/TestFunction.m
+++ b/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
diff --git a/lib/functions/TestFunctionGradient.m b/lib/functions/TestFunctionGradient.m
index 1ee153976e67724b6a9ba10b44c9447265dab0e4..8e9a94462cdc21ab3f6fc44c4136d1ee9610d629 100644
--- a/lib/functions/TestFunctionGradient.m
+++ b/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
diff --git a/lib/mesh/@Mesh/Mesh.m b/lib/mesh/@Mesh/Mesh.m
index 3305d23e1cfd6e3b4fe1b7657cc28d58c6e04be2..02333f174279cfbad7b491361d1ec0ab8b35d636 100644
--- a/lib/mesh/@Mesh/Mesh.m
+++ b/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.
@@ -75,8 +71,8 @@ classdef Mesh < handle
end
%% protected methods
- methods (Access='protected')
- updateData(obj)
+ methods (Access=protected)
+ updateData(obj, bisecData)
end
%% get methods for dependent data
diff --git a/lib/mesh/@Mesh/refineLocally.m b/lib/mesh/@Mesh/refineLocally.m
index 5fb042bcb45a33e2e60fe739397c8b466b1453ed..972ac8217456506865d55d3e6170f0c924c5ea2e 100644
--- a/lib/mesh/@Mesh/refineLocally.m
+++ b/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
diff --git a/lib/mesh/@Mesh/updateData.m b/lib/mesh/@Mesh/updateData.m
index 2756179f204255ecc375d63a604149ce116b6a3f..7cdad5624b01c306a8bbaed11847efa18633b69f 100644
--- a/lib/mesh/@Mesh/updateData.m
+++ b/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
-
-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
+end
\ No newline at end of file
diff --git a/lib/mesh/refinements/NVB.m b/lib/mesh/refinements/NVB.m
index 9a7384395d2e67e4ecdef01c707c486b42824c97..ae819ae629cedc8942fbd9a76156c5b0e6b4eb51 100644
--- a/lib/mesh/refinements/NVB.m
+++ b/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
diff --git a/lib/mesh/refinements/NVB1.m b/lib/mesh/refinements/NVB1.m
index 7af4bd8794401e163ea71a2fc0f7d92cb73aa7c1..59bf6bfcc9d3b698121a60e6455c1e1ea56c80eb 100644
--- a/lib/mesh/refinements/NVB1.m
+++ b/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;
diff --git a/lib/mesh/refinements/NVB3.m b/lib/mesh/refinements/NVB3.m
index a6a14a8ddfef1b66873485a1b8702b45451c98f4..7ab6e12a98f9b166f183cdd40776fac2983f2db5 100644
--- a/lib/mesh/refinements/NVB3.m
+++ b/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
\ No newline at end of file
+end
diff --git a/lib/mesh/refinements/NVB5.m b/lib/mesh/refinements/NVB5.m
index 15ff8de243e4aaa9bde81a5e8d95555b838e5123..aa8f4116cfb6cd8b31fd588a96df308e81417634 100644
--- a/lib/mesh/refinements/NVB5.m
+++ b/lib/mesh/refinements/NVB5.m
@@ -4,27 +4,30 @@
classdef NVB5 < NVB
%% methods
- methods (Access='public')
+ methods (Access=public)
function obj = NVB5(mesh)
obj = obj@NVB(mesh);
end
end
%% methods
- methods (Access='protected')
+ methods (Access=protected)
function markedEdges = markedElementsToEdges(obj, markedElements)
markedEdges = markedElementsToEdges@NVB(obj, markedElements);
end
- function bisecGroups = groupElements(obj, markedEdges, markedElements)
- bisecGroups = groupElements@NVB(obj, markedEdges);
-
+ function bisecData = groupElements(obj, markedEdges, markedElements)
% split elements with 3 marked edges into 'has interior node' and 'doesnt'
- noInteriorNode = false(1,obj.mesh.nElements);
- noInteriorNode(bisecGroups{4}.elementIdx) = true;
- noInteriorNode(markedElements) = false;
- bisecGroups{4} = Bisec123(find(noInteriorNode));
- bisecGroups{5} = Bisec5(markedElements);
+ hasInteriorNode = false(1,obj.mesh.nElements);
+ hasInteriorNode(markedElements) = true;
+
+ edge = reshape(markedEdges(obj.mesh.element2edges), size(obj.mesh.element2edges));
+ 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,:) & ~hasInteriorNode), ...
+ Bisec5, find(edge(1,:) & edge(2,:) & edge(3,:) & hasInteriorNode));
end
end
end
diff --git a/lib/mesh/refinements/NVBEdge.m b/lib/mesh/refinements/NVBEdge.m
index 735c87c5fb942be8a87f18fac05dc01fef179fc8..cf1734737f42152fd15b82354aab833d54914f34 100644
--- a/lib/mesh/refinements/NVBEdge.m
+++ b/lib/mesh/refinements/NVBEdge.m
@@ -4,14 +4,14 @@
classdef NVBEdge < NVB
%% methods
- methods (Access='public')
+ methods (Access=public)
function obj = NVBEdge(mesh)
obj = obj@NVB(mesh);
end
end
%% internal methods
- methods (Access='protected')
+ methods (Access=protected)
function markedEdges = markedElementsToEdges(obj, markedElements)
markedEdges = false(1,obj.mesh.nEdges);
markedEdges(markedElements) = true;
diff --git a/lib/mesh/refinements/RGB.m b/lib/mesh/refinements/RGB.m
index aa894349ed16c1da921cf57cef712a07a4f5e22c..cf8de10cb1234d927581efbeaf9fc71f4da78881 100644
--- a/lib/mesh/refinements/RGB.m
+++ b/lib/mesh/refinements/RGB.m
@@ -4,14 +4,14 @@
classdef RGB < NVB
%% methods
- methods (Access='public')
+ methods (Access=public)
function obj = RGB(mesh)
obj = obj@NVB(mesh);
end
end
%% methods
- methods (Access='protected')
+ methods (Access=protected)
function markedEdges = markedElementsToEdges(obj, markedElements)
% sort edges by length (while still retaining the orientation of the element)
mesh = obj.mesh;
@@ -23,9 +23,13 @@ classdef RGB < NVB
markedEdges(mesh.element2edges(:,markedElements)) = true;
end
- function bisecGroups = groupElements(obj, markedEdges, markedElements)
- bisecGroups = groupElements@NVB(obj, markedEdges, markedElements);
- bisecGroups{4} = BisecRed(bisecGroups{4}.elementIdx) ;
+ function bisecData = groupElements(obj, markedEdges, ~)
+ edge = reshape(markedEdges(obj.mesh.element2edges), size(obj.mesh.element2edges));
+ 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,:)), ...
+ BisecRed, find(edge(1,:) & edge(2,:) & edge(3,:)));
end
end
end
diff --git a/lib/mesh/refinements/bisection/AbstractBisection.m b/lib/mesh/refinements/bisection/AbstractBisection.m
index 7cd69ee80d5e3fa5f27abcb6c61b27afa6020453..11eed9282a363a6985bf1ebade8ffc7b7bb48a56 100644
--- a/lib/mesh/refinements/bisection/AbstractBisection.m
+++ b/lib/mesh/refinements/bisection/AbstractBisection.m
@@ -3,46 +3,21 @@
classdef AbstractBisection
%% properties
- properties (Abstract, GetAccess='public', SetAccess='protected')
+ properties (Abstract, GetAccess=public, SetAccess=protected)
innerNodes (3,:) double
nInnerEdges (1,1) double
nDescendants (1,1) double
end
- properties (Dependent)
- nInnerNodes
- nMembers
- end
-
- properties (GetAccess='public', SetAccess='protected')
- elementIdx (:,1) double
- end
-
%% methods
- methods (Access='public')
- function obj = AbstractBisection(elementIndices)
- obj.elementIdx = elementIndices;
- end
- end
-
- methods
- function n = get.nInnerNodes(obj)
- n = size(obj.innerNodes, Dim.Elements);
- end
-
- function n = get.nMembers(obj)
- n = length(obj.elementIdx);
- end
- end
-
- methods (Access='protected')
- function [T, F] = getTrueFalseArrays(obj)
- T = true(1, obj.nMembers);
- F = false(1, obj.nMembers);
+ methods (Access=protected)
+ function [T, F] = getTrueFalseArrays(obj, n)
+ T = true(1, n);
+ F = false(1, n);
end
end
- methods (Abstract, Access='public')
+ methods (Abstract, Access=public)
newElements = refineElement(obj, elements, edge2newNodes, innerNodes)
newEdges = createNewEdges(obj, elements, edge2newNodes, innerNodes)
newElement2Edges = refineEdgeConnectivity(obj, element2edges, left, right, innerEdges)
diff --git a/lib/mesh/refinements/bisection/Bisec1.m b/lib/mesh/refinements/bisection/Bisec1.m
index 0f371ff69ea4752d22e69c4e2bad8afcd55875f6..a600835b33e4637a558bc542d37d5211d0c12ef2 100644
--- a/lib/mesh/refinements/bisection/Bisec1.m
+++ b/lib/mesh/refinements/bisection/Bisec1.m
@@ -3,14 +3,14 @@
classdef Bisec1 < AbstractBisection
%% properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
innerNodes = []
nInnerEdges = 1
nDescendants = 2
end
%% methods
- methods (Access='public')
+ methods (Access=public)
function newElements = refineElement(~, oldNodes, edgeMidPts, ~)
newElements = [oldNodes(3,:), oldNodes(2,:); ...
oldNodes(1,:), oldNodes(3,:); ...
@@ -28,7 +28,7 @@ classdef Bisec1 < AbstractBisection
end
function newOrientation = refineEdgeOrientation(obj, old)
- [T, F] = getTrueFalseArrays(obj);
+ [T, F] = getTrueFalseArrays(obj, size(old,2));
newOrientation = [old(3,:), old(2,:); ...
old(1,:), T; ...
F, old(1,:)];
diff --git a/lib/mesh/refinements/bisection/Bisec12.m b/lib/mesh/refinements/bisection/Bisec12.m
index 01ddf9dadb1b548797b19a51c11d151487e54543..d39225ed6b17b70aa590aa89676e1968695fec09 100644
--- a/lib/mesh/refinements/bisection/Bisec12.m
+++ b/lib/mesh/refinements/bisection/Bisec12.m
@@ -3,14 +3,14 @@
classdef Bisec12 < AbstractBisection
%% properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
innerNodes = []
nInnerEdges = 2
nDescendants = 3
end
%% methods
- methods (Access='public')
+ methods (Access=public)
function newElements = refineElement(~, oldNodes, edgeMidPts, ~)
newElements = [oldNodes(3,:), edgeMidPts(1,:), oldNodes(3,:); ...
oldNodes(1,:), oldNodes(2,:), edgeMidPts(1,:); ...
@@ -30,7 +30,7 @@ classdef Bisec12 < AbstractBisection
end
function newOrientation = refineEdgeOrientation(obj, old)
- [T, F] = getTrueFalseArrays(obj);
+ [T, F] = getTrueFalseArrays(obj, size(old,2));
newOrientation = [old(3,:), old(1,:), T; ...
old(1,:), old(2,:), T; ...
F, F, old(2,:)];
diff --git a/lib/mesh/refinements/bisection/Bisec123.m b/lib/mesh/refinements/bisection/Bisec123.m
index 49ce4036efbc5c14b5adfc987481db7765972878..25de0c7419061efb60634cc8248490c637d3c325 100644
--- a/lib/mesh/refinements/bisection/Bisec123.m
+++ b/lib/mesh/refinements/bisection/Bisec123.m
@@ -3,14 +3,14 @@
classdef Bisec123 < AbstractBisection
%% properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
innerNodes = []
nInnerEdges = 3
nDescendants = 4
end
%% methods
- methods (Access='public')
+ methods (Access=public)
function newElements = refineElement(~, oldNodes, edgeMidPts, ~)
newElements = [edgeMidPts(1,:), oldNodes(1,:), edgeMidPts(1,:), oldNodes(3,:); ...
oldNodes(3,:), edgeMidPts(1,:), oldNodes(2,:), edgeMidPts(1,:); ...
@@ -29,7 +29,7 @@ classdef Bisec123 < AbstractBisection
end
function newOrientation = refineEdgeOrientation(obj, old)
- [T, F] = getTrueFalseArrays(obj);
+ [T, F] = getTrueFalseArrays(obj, size(old,2));
newOrientation = [F, old(1,:), old(1,:), T; ...
old(3,:), T, old(2,:), T; ...
F, old(3,:), F, old(2,:)];
diff --git a/lib/mesh/refinements/bisection/Bisec13.m b/lib/mesh/refinements/bisection/Bisec13.m
index 94418f917914195460d8506a8f431a27ad17a825..f67888fca295828d2881854e894d9a078b143237 100644
--- a/lib/mesh/refinements/bisection/Bisec13.m
+++ b/lib/mesh/refinements/bisection/Bisec13.m
@@ -3,14 +3,14 @@
classdef Bisec13 < AbstractBisection
%% properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
innerNodes = []
nInnerEdges = 2
nDescendants = 3
end
%% public methods
- methods (Access='public')
+ methods (Access=public)
function newElements = refineElement(~, oldNodes, edgeMidPts, ~)
newElements = [edgeMidPts(1,:), oldNodes(1,:), oldNodes(2,:); ...
oldNodes(3,:), edgeMidPts(1,:), oldNodes(3,:); ...
@@ -29,7 +29,7 @@ classdef Bisec13 < AbstractBisection
end
function newOrientation = refineEdgeOrientation(obj, old)
- [T, F] = getTrueFalseArrays(obj);
+ [T, F] = getTrueFalseArrays(obj, size(old,2));
newOrientation = [F, old(1,:), old(2,:); ...
old(3,:), T, T; ...
F, old(3,:), old(1,:)];
diff --git a/lib/mesh/refinements/bisection/Bisec5.m b/lib/mesh/refinements/bisection/Bisec5.m
index 412b348547146664b3c2628edc7cdec0d1dac3bb..2800f8ed9af6c46c94111aca22305a31a9a1db2c 100644
--- a/lib/mesh/refinements/bisection/Bisec5.m
+++ b/lib/mesh/refinements/bisection/Bisec5.m
@@ -3,7 +3,7 @@
classdef Bisec5 < AbstractBisection
%% properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
innerNodes = [1;1;2]/4;
nInnerEdges = 6
nDescendants = 6
@@ -11,7 +11,7 @@ classdef Bisec5 < AbstractBisection
end
%% methods
- methods (Access='public')
+ methods (Access=public)
function newElements = refineElement(~, oldNodes, edgeMidPts, innerNodes)
newElements = [edgeMidPts(3,:), oldNodes(3,:), oldNodes(1,:), edgeMidPts(1,:), edgeMidPts(2,:), edgeMidPts(1,:); ...
edgeMidPts(1,:), edgeMidPts(3,:), edgeMidPts(1,:), oldNodes(2,:), oldNodes(3,:), edgeMidPts(2,:); ...
@@ -30,7 +30,7 @@ classdef Bisec5 < AbstractBisection
end
function newOrientation = refineEdgeOrientation(obj, old)
- [T, F] = getTrueFalseArrays(obj);
+ [T, F] = getTrueFalseArrays(obj, size(old,2));
newOrientation = [T, old(3,:), old(1,:), old(1,:), old(2,:), T; ...
F, F, F, old(2,:), F, F; ...
T, T, old(3,:), F, T, T];
diff --git a/lib/mesh/refinements/bisection/BisecRed.m b/lib/mesh/refinements/bisection/BisecRed.m
index 0125fdb2dca4ac8f6c98e365cde7ad2823975123..7435ac39a021d1389010da12a3537ba017cb394a 100644
--- a/lib/mesh/refinements/bisection/BisecRed.m
+++ b/lib/mesh/refinements/bisection/BisecRed.m
@@ -3,14 +3,14 @@
classdef BisecRed < AbstractBisection
%% properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
innerNodes = []
nInnerEdges = 3
nDescendants = 4
end
%% methods
- methods (Access='public')
+ methods (Access=public)
function newElements = refineElement(~, oldNodes, edgeMidPts, ~)
newElements = [oldNodes(1,:), edgeMidPts(1,:), edgeMidPts(3,:), edgeMidPts(2,:); ...
edgeMidPts(1,:), oldNodes(2,:), edgeMidPts(2,:), edgeMidPts(3,:); ...
@@ -29,7 +29,7 @@ classdef BisecRed < AbstractBisection
end
function newOrientation = refineEdgeOrientation(obj, old)
- [T, F] = getTrueFalseArrays(obj);
+ [T, F] = getTrueFalseArrays(obj, size(old,2));
newOrientation = [old(1,:), old(1,:), F, T; ...
F, old(2,:), old(2,:), T; ...
old(3,:), F, old(3,:), T];
diff --git a/lib/mesh/utils/AffineTransformation.m b/lib/mesh/utils/AffineTransformation.m
index a4026e71aaf4c2868ebd614771217f810f9e5f81..2a93023c97983dc44dbce58db1e59e8d3f7577d5 100644
--- a/lib/mesh/utils/AffineTransformation.m
+++ b/lib/mesh/utils/AffineTransformation.m
@@ -9,7 +9,7 @@
classdef AffineTransformation < handle
%% properties
- properties (GetAccess = 'public', SetAccess = 'private')
+ properties (GetAccess=public, SetAccess=private)
DFinv (4,:) double
area (1,:) double
ds (1,:) double
diff --git a/lib/mesh/utils/BisectionEventData.m b/lib/mesh/utils/BisectionEventData.m
new file mode 100644
index 0000000000000000000000000000000000000000..35d39cc1f0b7bcf542813b53f28aa1dc88f21f03
--- /dev/null
+++ b/lib/mesh/utils/BisectionEventData.m
@@ -0,0 +1,75 @@
+% BisectionEventData (subclass of event.EventData) Encapsulates data used for
+% bisection.
+
+classdef (ConstructOnLoad) BisectionEventData < event.EventData
+ properties (GetAccess=public, SetAccess=protected)
+ nElements (1,1) double
+ bisectedEdges (:,1) logical
+ bisection (:,1) cell
+ nNonRefinedElements (1,1) double
+ nRefinedElements (:,1) double
+ nDescendants (:,1) double
+ nInnerNodes (:,1) double
+ nInnerEdges (:,1) double
+ end
+
+ properties (Access=protected)
+ nonRefinedElements
+ refinedElements
+ end
+
+ properties (Dependent)
+ nBisecMethods
+ end
+
+ methods (Access=public)
+ function obj = BisectionEventData(nElements, bisectedEdges, bisecMethods, refinedElements)
+ arguments
+ nElements
+ bisectedEdges
+ end
+ arguments (Repeating)
+ bisecMethods
+ refinedElements
+ end
+
+ obj.nElements = nElements;
+ obj.bisectedEdges = bisectedEdges;
+ obj.bisection = bisecMethods;
+ obj.refinedElements = refinedElements;
+
+ obj.nRefinedElements = cellfun(@(x) numel(x), refinedElements);
+ obj.nDescendants = cellfun(@(x) x.nDescendants, bisecMethods);
+ obj.nInnerNodes = cellfun(@(x) size(x.innerNodes,2), bisecMethods);
+ obj.nInnerEdges = cellfun(@(x) x.nInnerEdges, bisecMethods);
+
+ elements = true(nElements, 1);
+ elements(horzcat(obj.refinedElements{:})) = false;
+ obj.nonRefinedElements = find(elements);
+ obj.nNonRefinedElements = numel(obj.nonRefinedElements);
+ end
+
+ function children = getNChildrenPerElement(obj)
+ children = accumarray(horzcat(obj.refinedElements{:})', ...
+ repelem(obj.nDescendants, obj.nRefinedElements), ...
+ [obj.nElements, 1], [], 1);
+ end
+
+ % 0 = non-refined, 1,... = refined according to bisection{k}
+ function idx = getRefinedElementIdx(obj, k)
+ assert((0 <= k) && (k <= obj.nBisecMethods))
+ if k == 0
+ idx = obj.nonRefinedElements;
+ else
+ idx = obj.refinedElements{k};
+ end
+ end
+ end
+
+ methods
+ function val = get.nBisecMethods(obj)
+ val = numel(obj.bisection);
+ end
+ end
+end
+
diff --git a/lib/quadrature/barycentric/Barycentric.m b/lib/quadrature/barycentric/Barycentric.m
index 670652e44fa1c41a6151a57d9c5a3c79542497d2..9cc4f4a171b9cc7c6a1abf46c3789b352b2c885f 100644
--- a/lib/quadrature/barycentric/Barycentric.m
+++ b/lib/quadrature/barycentric/Barycentric.m
@@ -10,17 +10,17 @@
classdef Barycentric
%% properties
- properties (GetAccess = 'public', SetAccess = 'protected')
+ properties (GetAccess=public, SetAccess=protected)
coordinates
nNodes
end
- properties (Abstract, GetAccess = 'public', SetAccess = 'protected')
+ properties (Abstract, GetAccess=public, SetAccess=protected)
dim
end
%% methods
- methods (Access = 'public')
+ methods (Access=public)
function obj = Barycentric(coordinates)
arguments
coordinates {Barycentric.verify}
diff --git a/lib/quadrature/barycentric/Barycentric1D.m b/lib/quadrature/barycentric/Barycentric1D.m
index 77379fe1de0c2867c1c569ee4ec8681ba841b100..d2cadd38404643669c95f52c978b2a4692d2ac8c 100644
--- a/lib/quadrature/barycentric/Barycentric1D.m
+++ b/lib/quadrature/barycentric/Barycentric1D.m
@@ -2,7 +2,7 @@
% coordinates on edges.
classdef Barycentric1D < Barycentric
- properties (GetAccess = 'public', SetAccess = 'protected')
+ properties (GetAccess=public, SetAccess=protected)
dim = 1;
end
-end
\ No newline at end of file
+end
diff --git a/lib/quadrature/barycentric/Barycentric2D.m b/lib/quadrature/barycentric/Barycentric2D.m
index b47e38d9fcbcfa0e00bffe8a30b523e9d78d7d30..9fa6444d00cb92448790a193774670a0059bb212 100644
--- a/lib/quadrature/barycentric/Barycentric2D.m
+++ b/lib/quadrature/barycentric/Barycentric2D.m
@@ -5,7 +5,7 @@
% Barycentric1D on edge given by edgeInd.
classdef Barycentric2D < Barycentric
- properties (GetAccess = 'public', SetAccess = 'protected')
+ properties (GetAccess=public, SetAccess=protected)
dim = 2;
end
@@ -24,4 +24,4 @@ classdef Barycentric2D < Barycentric
bary2D = Barycentric2D(newNodes);
end
end
-end
\ No newline at end of file
+end
diff --git a/lib/quadrature/quadratureRules/QuadratureRule.m b/lib/quadrature/quadratureRules/QuadratureRule.m
index 3fafec1c1a2ad26265b29947cc272e3d875e4305..38d2340340b1e091e0ce185a3ff0d45e65f630ad 100644
--- a/lib/quadrature/quadratureRules/QuadratureRule.m
+++ b/lib/quadrature/quadratureRules/QuadratureRule.m
@@ -14,7 +14,7 @@
classdef QuadratureRule
%% properties
- properties (GetAccess = 'public', SetAccess = 'protected')
+ properties (GetAccess=public, SetAccess=protected)
bary
weights
nNodes
@@ -22,7 +22,7 @@ classdef QuadratureRule
end
%% methods
- methods (Access = 'public')
+ methods (Access=public)
function obj = QuadratureRule(bary, weights)
arguments
bary (1,1) Barycentric
diff --git a/lib/spaces/@FeSpace/FeSpace.m b/lib/spaces/@FeSpace/FeSpace.m
index 1b8b72280172ad67311df65002b4d9a7c9eb628e..d0fbc2da0753f1f10c823febdc99020bd0488930 100644
--- a/lib/spaces/@FeSpace/FeSpace.m
+++ b/lib/spaces/@FeSpace/FeSpace.m
@@ -4,24 +4,24 @@
classdef FeSpace < handle
%% properties
- properties (GetAccess = 'public', SetAccess = 'protected')
+ properties (GetAccess=public, SetAccess=protected)
mesh
finiteElement
dirichlet
end
- properties (Access = 'private')
+ properties (Access=private)
dofs
freeDofs
end
%% events
- properties (Access = 'private')
+ properties (Access=private)
listenerHandle
end
%% methods
- methods (Access = 'public')
+ methods (Access=public)
function obj = FeSpace(mesh, fe, opt)
% Construct FeSpace from given mesh and finite element.
%
@@ -39,7 +39,7 @@ classdef FeSpace < handle
end
% register fespace as an observer of mesh
- obj.listenerHandle = mesh.listener('HasChanged', @obj.handleChangedMesh);
+ obj.listenerHandle = mesh.listener('JustRefined', @obj.handleChangedMesh);
% initialize lazily evaluated properties
obj.dofs = [];
@@ -79,12 +79,12 @@ classdef FeSpace < handle
end
end
- methods (Access = 'protected')
+ methods (Access=protected)
dofs = assembleDofs(obj)
freeDofs = assembleFreeDofs(obj)
end
- methods (Access = 'private')
+ methods (Access=private)
handleChangedMesh(obj, src, event)
end
end
diff --git a/lib/spaces/finiteElements/BernsteinBezierPoly.m b/lib/spaces/finiteElements/BernsteinBezierPoly.m
index d3bc3e581d46b753e2230a2307857227bcb0c082..e0cde89867db6cfab1442fa41980d243a42f0dcc 100644
--- a/lib/spaces/finiteElements/BernsteinBezierPoly.m
+++ b/lib/spaces/finiteElements/BernsteinBezierPoly.m
@@ -7,14 +7,14 @@
classdef BernsteinBezierPoly
%% properties
- properties (Access='protected')
+ properties (Access=protected)
constants
exponents
nFun
end
%% methods
- methods (Access='public')
+ methods (Access=public)
function obj = setExponents(obj, exponents)
arguments
obj
diff --git a/lib/spaces/finiteElements/FiniteElement.m b/lib/spaces/finiteElements/FiniteElement.m
index f8e313ba851aab13150af276fc358f845e34d92b..8ae7891630d9672b4737d681591149b0e751431a 100644
--- a/lib/spaces/finiteElements/FiniteElement.m
+++ b/lib/spaces/finiteElements/FiniteElement.m
@@ -3,7 +3,7 @@
% the dofs are connected.
classdef FiniteElement
- methods (Abstract, Access = 'public')
+ methods (Abstract, Access=public)
% evalShapeFunctions Evaluate local shapefunctions at barycentric
% coordinates
evalShapeFunctions(obj, bary)
diff --git a/lib/spaces/finiteElements/HigherOrderH1Fe.m b/lib/spaces/finiteElements/HigherOrderH1Fe.m
index cfaa48fdd57d5c1f34ce1c60b54d303a4fceefd2..2ce914d9472d45100fa6d3f5e871fc03eddb284d 100644
--- a/lib/spaces/finiteElements/HigherOrderH1Fe.m
+++ b/lib/spaces/finiteElements/HigherOrderH1Fe.m
@@ -4,7 +4,7 @@
classdef HigherOrderH1Fe < NodalFiniteElement & BernsteinBezierPoly
%% implementation of abstract superclass properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
order
nodeDofLocations
edgeDofLocations
@@ -12,7 +12,7 @@ classdef HigherOrderH1Fe < NodalFiniteElement & BernsteinBezierPoly
end
%% implementation of abstract superclass methods
- methods (Access = 'public')
+ methods (Access=public)
function obj = HigherOrderH1Fe(order)
arguments
order (1,1) double {mustBePositive}
diff --git a/lib/spaces/finiteElements/HigherOrderL2Fe.m b/lib/spaces/finiteElements/HigherOrderL2Fe.m
index 9db6258291fc916f7e4addf71d822b1f5e300bc3..c8dd2e049a36d73138c869f22fad840a7457d2de 100644
--- a/lib/spaces/finiteElements/HigherOrderL2Fe.m
+++ b/lib/spaces/finiteElements/HigherOrderL2Fe.m
@@ -4,7 +4,7 @@
classdef HigherOrderL2Fe < NodalFiniteElement & BernsteinBezierPoly
%% implementation of abstract superclass properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
order
nodeDofLocations
edgeDofLocations
@@ -12,7 +12,7 @@ classdef HigherOrderL2Fe < NodalFiniteElement & BernsteinBezierPoly
end
%% implementation of abstract superclass methods
- methods (Access = 'public')
+ methods (Access=public)
%% constructor
function obj = HigherOrderL2Fe(order)
arguments
diff --git a/lib/spaces/finiteElements/LowestOrderH1Fe.m b/lib/spaces/finiteElements/LowestOrderH1Fe.m
index 8cde03cdb4ae2c70e09fa4ed565dddc5cc29309a..7bfe66c2d320c66800a5aa49da3cf7c1c77a3823 100644
--- a/lib/spaces/finiteElements/LowestOrderH1Fe.m
+++ b/lib/spaces/finiteElements/LowestOrderH1Fe.m
@@ -4,7 +4,7 @@
classdef LowestOrderH1Fe < NodalFiniteElement
%% implementation of abstract superclass properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
order
nodeDofLocations
edgeDofLocations
@@ -12,7 +12,7 @@ classdef LowestOrderH1Fe < NodalFiniteElement
end
%% implementation of abstract superclass methods
- methods (Access = 'public')
+ methods (Access=public)
function obj = LowestOrderH1Fe()
obj.order = 1;
obj.nodeDofLocations = 1;
diff --git a/lib/spaces/finiteElements/LowestOrderL2Fe.m b/lib/spaces/finiteElements/LowestOrderL2Fe.m
index 7b5002c7a40cb92a00510f64f0d37ba93c12e1df..8c11644e983077cfb1811e65de85f3518f94bf08 100644
--- a/lib/spaces/finiteElements/LowestOrderL2Fe.m
+++ b/lib/spaces/finiteElements/LowestOrderL2Fe.m
@@ -4,7 +4,7 @@
classdef LowestOrderL2Fe < NodalFiniteElement
%% implementation of abstract superclass properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
order
nodeDofLocations
edgeDofLocations
@@ -12,7 +12,7 @@ classdef LowestOrderL2Fe < NodalFiniteElement
end
%% implementation of abstract superclass methods
- methods (Access = 'public')
+ methods (Access=public)
function obj = LowestOrderL2Fe()
obj.order = 0;
obj.nodeDofLocations = [];
diff --git a/lib/spaces/finiteElements/NodalFiniteElement.m b/lib/spaces/finiteElements/NodalFiniteElement.m
index 58da69af0cc6e7bc8484dd2a90959244d91f3488..9766fd53e54ffed82cd6687d021c3a2bd6f52be7 100644
--- a/lib/spaces/finiteElements/NodalFiniteElement.m
+++ b/lib/spaces/finiteElements/NodalFiniteElement.m
@@ -5,7 +5,7 @@
classdef NodalFiniteElement < FiniteElement
%% properties
- properties (Abstract, GetAccess='public', SetAccess='protected')
+ properties (Abstract, GetAccess=public, SetAccess=protected)
nodeDofLocations
edgeDofLocations
innerDofLocations
diff --git a/lib/spaces/interpolation/FeProlongation.m b/lib/spaces/interpolation/FeProlongation.m
new file mode 100644
index 0000000000000000000000000000000000000000..4dc5d6b97ac4e223403a4794c378ad118f6e3a8a
--- /dev/null
+++ b/lib/spaces/interpolation/FeProlongation.m
@@ -0,0 +1,134 @@
+% FeProlongation (subclass of Prolongation) Provides prolongation operator for
+% general FeFunction to a refined mesh.
+%
+% P = FeProlongation(fes) returns a handle to the prolongation object
+% associated to the finite element space fes. The prolongation matrix
+% P.matrix is set automatically at mesh refinement.
+%
+% prolongate(P, u) returns the prolongated data of FeFunction u.
+
+classdef FeProlongation < Prolongation
+ %% properties
+ properties (Access=private)
+ postRefineListener
+ end
+
+ %% methods
+ methods (Access=public)
+ function obj = FeProlongation(fes)
+ obj = obj@Prolongation(fes);
+ obj.postRefineListener = fes.mesh.listener('RefineCompleted', @obj.connectDofs);
+
+ assert(isa(obj.fes.finiteElement, 'NodalFiniteElement'), ...
+ 'Prolongation defined only for nodal finite elements.')
+ end
+ end
+
+ methods (Access=protected)
+ function setupMatrix(obj, ~, data)
+ % general idea: compute *all* dofs for each new element and store
+ % them consecutively
+ % those are connected to the actual new dofs when their numbering
+ % per element is available (-> obj.connectDofs)
+
+ % de-allocate for memory-efficiency
+ obj.matrix = [];
+
+ % get dof information
+ oldDofs = getDofs(obj.fes);
+ fe = obj.fes.finiteElement;
+ dofLocations = getDofLocations(fe);
+ nLocalDofs = dofLocations.nNodes;
+
+ % pre-allocate
+ element2newDof = cumsum([1; getNChildrenPerElement(data)*nLocalDofs]);
+ n = data.nNonRefinedElements*nLocalDofs + sum(data.nRefinedElements)*nLocalDofs^2;
+ [I, J, V] = deal(zeros(n, 1));
+
+ % for non-refined elements, just transfer the dofs
+ n = data.nNonRefinedElements*nLocalDofs;
+ if n ~= 0
+ idx = (1:n)';
+ nonRefinedElems = getRefinedElementIdx(data, 0);
+ I(idx) = getConsecutiveIndices(element2newDof(nonRefinedElems), nLocalDofs);
+ J(idx) = reshape(oldDofs.element2Dofs(:,nonRefinedElems), [], 1);
+ V(idx) = 1;
+ end
+ ptr = n;
+
+ % handle refined elements (nodal interpolation on new elements)
+ localMatrix = squeeze(evalShapeFunctions(fe, dofLocations));
+ for k = find(data.nRefinedElements ~= 0)'
+ unitTriangle = getBisectedUnitTriangle(class(data.bisection{k}));
+ newDofLocations = getNewLocationsElementwise(unitTriangle, dofLocations);
+ newDofWeights = divideElementwise(...
+ reshape(evalShapeFunctions(fe, newDofLocations), nLocalDofs, nLocalDofs, []), localMatrix);
+
+ nNewLocs = newDofLocations.nNodes;
+ elems = getRefinedElementIdx(data, k);
+ n = data.nRefinedElements(k) * numel(newDofWeights);
+ idx = ptr + (1:n)';
+ I(idx) = repelem(getConsecutiveIndices(element2newDof(elems), nNewLocs), nLocalDofs);
+ J(idx) = reshape(repelem(oldDofs.element2Dofs(:,elems), 1, nNewLocs), [], 1);
+ V(idx) = reshape(repmat(newDofWeights, 1, data.nRefinedElements(k)), [], 1);
+ ptr = ptr + n;
+ end
+
+ nNewDofs = element2newDof(end)-1;
+ obj.matrix = sparse(I, J, V, nNewDofs, oldDofs.nDofs);
+ end
+
+ function connectDofs(obj, ~, ~)
+ dofs = getDofs(obj.fes);
+ nOutput = size(obj.matrix, 1);
+ idx = zeros(dofs.nDofs, 1);
+
+ idx(dofs.element2Dofs(:)) = 1:nOutput;
+ obj.matrix = obj.matrix(idx,:);
+ end
+ end
+end
+
+%% local functions
+function mesh = getBisectedUnitTriangle(bisecMethod)
+ mesh = Mesh.loadFromGeometry('unittriangle');
+
+ % bisect unit triangle according to given bisection rule
+ % TODO: this should be handled by bisection method itself?
+ switch bisecMethod
+ case 'Bisec1'
+ mesh.refineLocally(3, 'NVBEdge');
+ case 'Bisec12'
+ mesh.refineLocally([2,3], 'NVBEdge');
+ case 'Bisec13'
+ mesh.refineLocally([1,3], 'NVBEdge');
+ case 'Bisec123'
+ mesh.refineUniform(1, 'NVB');
+ case 'Bisec5'
+ mesh.refineUniform(1, 'NVB5');
+ case 'BisecRed'
+ mesh.refineUniform(1, 'RGB');
+ otherwise
+ error(['Bisection method ', bisecMethod, ' not related to refinement of unit triangle.'])
+ end
+end
+
+function elementwiseLocations = getNewLocationsElementwise(mesh, locations)
+ % collect dofs such that all locations of one child appear before those of other children
+ % -> permute dimensions: 2 = children, 3 = barycentric coordinates
+ elementwiseLocations = permute(elementwiseCoordinates(mesh, locations), [1,3,2]);
+ elementwiseLocations = reshape(elementwiseLocations, 2, [], 1);
+ elementwiseLocations = Barycentric2D(...
+ cat(Dim.Vector, elementwiseLocations, 1-sum(elementwiseLocations, Dim.Vector)));
+end
+
+function C = divideElementwise(A, B)
+ % TODO: from 2022a on, this can be replaced by pagemrdivide(A, B)
+ n = size(A,1);
+ C = B' \ reshape(permute(A, [2,1,3]), n, []);
+ C = reshape(permute(reshape(C, n, n, []), [2,1,3]), n, []);
+end
+
+function idx = getConsecutiveIndices(parents, nIndices)
+ idx = reshape(((0:(nIndices-1)) + parents)', [], 1);
+end
diff --git a/lib/spaces/interpolation/LinearFeProlongation.m b/lib/spaces/interpolation/LinearFeProlongation.m
deleted file mode 100644
index b28782363a7472016fa0b430b32aed57f0b24f1a..0000000000000000000000000000000000000000
--- a/lib/spaces/interpolation/LinearFeProlongation.m
+++ /dev/null
@@ -1,43 +0,0 @@
-% LinearFeProlongation (subclass of Prolongation) Provides piecewise linear
-% interpolant of FeFunction on a refined mesh.
-%
-% Pu = LinearFeProlongation(u) returns the piecewise linear interpolation of u
-% that can be prolongated to a finer mesh.
-%
-% Pu.interpolatedData is the data of the prolongation and is set automatically
-% at mesh refinement.
-
-classdef LinearFeProlongation < Prolongation
- %% properties
- properties (GetAccess='public', SetAccess='protected')
- interpolatedData
- end
-
- properties (Access='private')
- interpolant
- end
-
- %% methods
- methods (Access='public')
- function obj = LinearFeProlongation(feFunction)
- obj = obj@Prolongation(feFunction);
- obj.interpolant = [];
- end
- end
-
- methods
- function val = get.interpolatedData(obj)
- if ~isempty(obj.interpolant)
- obj.interpolatedData = nodalInterpolation(obj.interpolant, obj.u.fes);
- obj.interpolant = [];
- end
- val = obj.interpolatedData;
- end
- end
-
- methods (Access='protected')
- function interpolate(obj, ~, ~)
- obj.interpolant = LinearInterpolation(obj.u);
- end
- end
-end
diff --git a/lib/spaces/interpolation/LinearInterpolation.m b/lib/spaces/interpolation/LinearInterpolation.m
index e5482a882736d1c817fa4ac8b261229082d10286..d5d1e7db815f156d168e6d17cb4d3fe2a0844ab2 100644
--- a/lib/spaces/interpolation/LinearInterpolation.m
+++ b/lib/spaces/interpolation/LinearInterpolation.m
@@ -5,17 +5,17 @@
classdef LinearInterpolation < Evaluable
%% properties
- properties (GetAccess='public', SetAccess='protected')
+ properties (GetAccess=public, SetAccess=protected)
mesh
u
end
- properties (Access='private')
+ properties (Access=private)
interpolant
end
%% methods
- methods (Access='public')
+ methods (Access=public)
function obj = LinearInterpolation(u)
arguments
u Evaluable
diff --git a/lib/spaces/interpolation/LoFeProlongation.m b/lib/spaces/interpolation/LoFeProlongation.m
new file mode 100644
index 0000000000000000000000000000000000000000..2c2603e8cf42a7371dabb010b93874d389047f30
--- /dev/null
+++ b/lib/spaces/interpolation/LoFeProlongation.m
@@ -0,0 +1,91 @@
+% LoFeProlongation (subclass of Prolongation) Prolongate lowest order L^2/H^1
+% conforming finite element function to refined mesh.
+%
+% P = LoFeProlongation(fes) returns a handle to the prolongation object
+% associated to the finite element space fes. The prolongation matrix
+% P.matrix is set automatically at mesh refinement.
+%
+% prolongate(P, u) returns the prolongated data of FeFunction u.
+
+classdef LoFeProlongation < Prolongation
+ %% properties
+ properties (Access=protected)
+ feType
+ end
+ %% methods
+ methods (Access=public)
+ function obj = LoFeProlongation(fes)
+ obj = obj@Prolongation(fes);
+ switch class(fes.finiteElement)
+ case 'LowestOrderL2Fe'
+ obj.feType = 'L2';
+ case 'LowestOrderH1Fe'
+ obj.feType = 'H1';
+ otherwise
+ eid = 'LoFeProlongation:wrongFeType';
+ msg = 'LoFeProlongation needs a lowest order L2 or H1 finite element space.';
+ throwAsCaller(MException(eid, msg));
+ end
+ end
+ end
+
+ methods (Access=protected)
+ function setupMatrix(obj, mesh, data)
+ if isequal(obj.feType, 'L2')
+ setupMatrixL2(obj, mesh, data);
+ elseif isequal(obj.feType, 'H1')
+ setupMatrixH1(obj, mesh, data);
+ end
+ end
+
+ function setupMatrixL2(obj, mesh, data)
+ % use that for lowest order L2 elements the dofs correspond to
+ % elements and indices of new elements is known
+ nChildren = getNChildrenPerElement(data);
+ nNewElements = sum(nChildren);
+
+ obj.matrix = sparse(1:nNewElements, repelem(1:mesh.nElements, nChildren), ...
+ ones(nNewElements,1), nNewElements, mesh.nElements);
+ end
+
+ function setupMatrixH1(obj, mesh, data)
+ % use that for lowest order H1 elements the dofs correspond to
+ % coordinates and new coordinates reside on edges or on inner nodes
+ dofs = getDofs(obj.fes);
+ nEntries = mesh.nCoordinates + 2*nnz(data.bisectedEdges) + ...
+ 3*sum(data.nInnerNodes.*data.nRefinedElements);
+ [I, J, V] = deal(zeros(nEntries, 1));
+
+ % node dofs stay the same
+ idx = 1:mesh.nCoordinates;
+ I(idx) = idx';
+ J(idx) = idx';
+ V(idx) = ones(idx(end), 1);
+ dofNr = mesh.nCoordinates;
+
+ % edge dofs are mean of edge dofs
+ n = nnz(data.bisectedEdges);
+ idx = idx(end) + (1:2*n);
+ I(idx) = repelem(dofNr + (1:n)', 2);
+ J(idx) = reshape(dofs.edge2Dofs(:,data.bisectedEdges), [], 1);
+ V(idx) = ones(2*n, 1)/2;
+ dofNr = idx(end) + n;
+
+ % inner dofs are weighted sum of element dofs
+ idxEnd = idx(end);
+ for k = find(data.nRefinedElements)'
+ n = data.nRefinedElements(k)*data.nInnerNodes(k);
+ idx = idxEnd + (1:3*n);
+ I(idx) = repelem(dofNr + (1:n)', 3);
+ J(idx) = reshape(repelem(dofs.element2Dofs, [1;1;1], data.nInnerNodes(k)), [], 1);
+ V(idx) = reshape(repmat(data.bisection{k}.innerNodes, 1, data.nInnerNodes(k)), [], 1);
+ dofNr = dofNr + n;
+ idxEnd = idxEnd + 3*n;
+ end
+
+ nNewDofs = mesh.nCoordinates + nnz(data.bisectedEdges) + ...
+ sum(data.nInnerNodes.*data.nRefinedElements);
+ obj.matrix = sparse(I, J, V, nNewDofs, dofs.nDofs);
+ end
+ end
+end
diff --git a/lib/spaces/interpolation/LoH1Prolongation.m b/lib/spaces/interpolation/LoH1Prolongation.m
deleted file mode 100644
index 10e505accb35679b03c72e65bd825074140b8c70..0000000000000000000000000000000000000000
--- a/lib/spaces/interpolation/LoH1Prolongation.m
+++ /dev/null
@@ -1,41 +0,0 @@
-% LoH1Prolongation (subclass of Prolongation) Prolongate lowest order H^1
-% conforming finite element function to refined mesh.
-%
-% Pu = LoH1Prolongation(u) returns a handle to the prolongation object.
-%
-% Pu.interpolatedData is the data of the prolongation and is set automatically
-% at mesh refinement.
-
-classdef LoH1Prolongation < Prolongation
- %% properties
- properties (GetAccess='public', SetAccess='protected')
- interpolatedData
- end
-
- %% methods
- methods (Access='public')
- function obj = LoH1Prolongation(feFunction)
- assert(isa(feFunction.fes.finiteElement, 'LowestOrderH1Fe'), ...
- 'Prolongation only implemented for lowest order H1 conforming finite elements.')
- obj = obj@Prolongation(feFunction);
- end
- end
-
- methods (Access='protected')
- function interpolate(obj, src, ~)
- % use that for lowest order H1 elements the dofs correspond to
- % coordinates and new coordinates reside on edges or on inner nodes
- uInnerNodes = cell(numel(src.bisecGroups), 1);
- for k = 1:numel(src.bisecGroups)
- if ~isempty(src.bisecGroups{k}.innerNodes)
- uInnerNodes{k} = eval(obj.u, src.bisecGroups{k}.innerNodes, ...
- src.bisecGroups{k}.elementIdx);
- end
- end
- midpoint = Barycentric1D([1;1]/2);
- obj.interpolatedData = [obj.u.data, ...
- evalEdge(obj.u, midpoint, src.markedEdges), ...
- horzcat(uInnerNodes{:})];
- end
- end
-end
diff --git a/lib/spaces/interpolation/LoL2Prolongation.m b/lib/spaces/interpolation/LoL2Prolongation.m
deleted file mode 100644
index a3928cf7968921e404147eb1248bf0e8411914c6..0000000000000000000000000000000000000000
--- a/lib/spaces/interpolation/LoL2Prolongation.m
+++ /dev/null
@@ -1,35 +0,0 @@
-% LoL2Prolongation (subclass of Prolongation) Prolongate lowest order L^2 finite
-% element function to refined mesh.
-%
-% Pu = LoL2Prolongation(u) returns a handle to the prolongation object.
-%
-% Pu.interpolatedData is the data of the prolongation and is set automatically
-% at mesh refinement.
-
-classdef LoL2Prolongation < Prolongation
- %% properties
- properties (GetAccess='public', SetAccess='protected')
- interpolatedData
- end
-
- %% methods
- methods (Access='public')
- function obj = LoL2Prolongation(feFunction)
- assert(isa(feFunction.fes.finiteElement, 'LowestOrderL2Fe'), ...
- 'Prolongation only implemented for lowest order L2 finite elements.')
- obj = obj@Prolongation(feFunction);
- end
- end
-
- methods (Access='protected')
- function interpolate(obj, src, ~)
- % use that for lowest order L2 elements the dofs correspond to elements
- % and indices of new elements is known
- nChildElems = ones(src.nElements, 1);
- for k = 1:length(src.bisecGroups)
- nChildElems(src.bisecGroups{k}.elementIdx) = src.bisecGroups{k}.nDescendants;
- end
- obj.interpolatedData = repelem(obj.u.data, nChildElems);
- end
- end
-end
diff --git a/lib/spaces/interpolation/Prolongation.m b/lib/spaces/interpolation/Prolongation.m
index 709b535fe0a0026c662318e91fbd2ffa1f5b0ec2..97adb20868910b445d6a1884181623ff17abeab0 100644
--- a/lib/spaces/interpolation/Prolongation.m
+++ b/lib/spaces/interpolation/Prolongation.m
@@ -3,26 +3,46 @@
classdef Prolongation < handle
%% properties
- properties (Abstract, GetAccess='public', SetAccess='protected')
- interpolatedData
+ properties (GetAccess=public, SetAccess=protected)
+ matrix {mustBeSparse}
end
- properties (Access='protected')
- u FeFunction
+ properties (Access=protected)
+ fes FeSpace
listenerHandle
end
%% methods
- methods (Access='public')
- function obj = Prolongation(feFunction)
- obj.u = feFunction;
- mesh = feFunction.fes.mesh;
- obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.interpolate);
- obj.interpolatedData = [];
+ methods (Access=public)
+ function obj = Prolongation(fes)
+ obj.fes = fes;
+ mesh = fes.mesh;
+ obj.matrix = [];
+ obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.setupMatrix);
+ end
+
+ function newFeData = prolongate(obj, u)
+ arguments
+ obj
+ u FeFunction
+ end
+
+ if isempty(obj.matrix)
+ error('Prolongation matrix not set. Are you sure there was mesh refinement?')
+ end
+ newFeData = obj.matrix * u.data';
end
end
- methods (Abstract, Access='protected')
- interpolate(obj, src, event)
+ methods (Abstract, Access=protected)
+ setupMatrix(obj, src, event)
+ end
+end
+
+function mustBeSparse(A)
+ if ~(isempty(A) || issparse(A))
+ eidType = 'Prolongation:matrixNotSparse';
+ msgType = 'Property ''matrix'' must be a sparse matrix.';
+ throwAsCaller(MException(eidType,msgType))
end
end
diff --git a/tests/TestFiniteElementFunctions.m b/tests/TestFiniteElementFunctions.m
index 9ae9f7365dced10872448045243e930c00f62494..f208227e5759ddeef015997f8a24a3eea1b81e43 100644
--- a/tests/TestFiniteElementFunctions.m
+++ b/tests/TestFiniteElementFunctions.m
@@ -66,7 +66,7 @@ methods (Test)
testCase.verifyEqual(nnz(testCase.u.data == scalar), nFreeDofs)
end
- function linearInterpolationPreservesConstants(testCase)
+ function prolongationPreservesConstants(testCase)
mesh = Mesh.loadFromGeometry('unitsquare');
fes = FeSpace(mesh, testCase.u.fes.finiteElement);
v = FeFunction(fes);
@@ -75,20 +75,42 @@ methods (Test)
integral = sum(integrateElement(v, qr), Dim.Elements);
testCase.verifyEqual(integral, pi, 'AbsTol', sqrt(eps))
end
+
+ function generalProlongationPreservesBasisFunctions(testCase)
+ mesh = Mesh.loadFromGeometry('unitsquare');
+ fes = FeSpace(mesh, testCase.u.fes.finiteElement);
+ qr = QuadratureRule.ofOrder(max(fes.finiteElement.order,1));
+ v = FeFunction(fes);
+ P = FeProlongation(fes);
+ for k = 1:min(5, size(getDofs(fes).element2Dofs, 1))
+ testCase.setToElementwiseBasisFunction(v, k);
+ before = sum(integrateElement(v, qr));
+ mesh.refineLocally(1, 'NVB');
+ v.setData(prolongate(P, v));
+ after = sum(integrateElement(v, qr));
+ testCase.verifyEqual(before, after, 'AbsTol', sqrt(eps))
+ end
+ end
end
-methods (Access='private')
+methods (Access=private)
function refineAndInterpolate(~, mesh, fes, v, val)
v.setData(val);
- if isa(fes.finiteElement, 'LowestOrderH1Fe')
- oldVals = LoH1Prolongation(v);
- elseif isa(fes.finiteElement, 'LowestOrderL2Fe')
- oldVals = LoL2Prolongation(v);
+ if isa(fes.finiteElement, 'LowestOrderH1Fe') ...
+ || isa(fes.finiteElement, 'LowestOrderL2Fe')
+ P = LoFeProlongation(fes);
else
- oldVals = LinearFeProlongation(v);
+ P = FeProlongation(fes);
end
mesh.refineUniform();
- v.setData(oldVals.interpolatedData);
+ v.setData(prolongate(P, v));
+ end
+
+ function setToElementwiseBasisFunction(~, v, k)
+ v.setData(0);
+ data = v.data;
+ data(getDofs(v.fes).element2Dofs(k,:)) = 1;
+ v.setData(data);
end
end
diff --git a/tests/TestRefinement.m b/tests/TestRefinement.m
index 10d5b1c687b9e0c97ebe0977fcd8aa048f3f363f..3ce2b4dd249987a775612ee2db1197f7d293fe0c 100644
--- a/tests/TestRefinement.m
+++ b/tests/TestRefinement.m
@@ -30,7 +30,7 @@ methods (Test)
end
end
-methods (Access='protected')
+methods (Access=protected)
function verifyEdgeConnectivityIntegrity(testCase, mesh)
[edgesRef, e2eRef, ~, bndRef] = Mesh.computeEdgeInformation( ...
mesh.elements, cellfun(@(bnd) mesh.edges(:,bnd), mesh.boundaries, 'UniformOutput', false));
@@ -53,4 +53,4 @@ methods (Access='protected')
end
end
-end
\ No newline at end of file
+end
diff --git a/tests/TestVectorProduct.m b/tests/TestVectorProduct.m
index b0d55aad91f720be538b33dab725848569b98fb6..ba515e96636bc3a61d4f96bd781d91b9d8940572 100644
--- a/tests/TestVectorProduct.m
+++ b/tests/TestVectorProduct.m
@@ -43,7 +43,7 @@ methods (Test, ParameterCombination='exhaustive')
end
end
-methods (Access='protected')
+methods (Access=protected)
function [Aext, dimA] = extendAlongDim(~, A, dim, n)
dimA = size(A);
repetitions = ones(1, dim);
@@ -52,4 +52,4 @@ methods (Access='protected')
end
end
-end
\ No newline at end of file
+end