From a55003303c92fdb817340834c25de0252ab55067 Mon Sep 17 00:00:00 2001 From: Julian Streitberger <julian.streitberger@tuwien.ac.at> Date: Tue, 28 Jan 2025 17:15:58 +0100 Subject: [PATCH] Work in progress --- experiments/jumpingDiffusion.m | 80 ++++++++++++++++++++++++++++++ lib/functions/piecewiseEvaluable.m | 44 +++++++++++----- 2 files changed, 111 insertions(+), 13 deletions(-) create mode 100644 experiments/jumpingDiffusion.m diff --git a/experiments/jumpingDiffusion.m b/experiments/jumpingDiffusion.m new file mode 100644 index 0000000..73c9c28 --- /dev/null +++ b/experiments/jumpingDiffusion.m @@ -0,0 +1,80 @@ +mesh = Mesh.loadFromGeometry('unitsquare'); +mesh.refineUniform(1); +fes = FeSpace(mesh, HigherOrderH1Fe(1), 'dirichlet', ':'); + +%% set problem data +blf = BilinearForm(); + +ncFes = FeSpace(mesh, LowestOrderL2Fe()); +a = FeFunction(ncFes); +trafo = getAffineTransformation(mesh); +n = trafo.n; +f(mesh.coordinates, n) +%b = MeshFunction(mesh, @(x) x) +%b = PiecewiseEvaluable(mesh, @(x, nu) f(x, nu), n); +%a.setData(nodalInterpolation(MeshFunction(mesh, @(x) f(x, n)), ncFes)); +%blf.a = a; +%qr = QuadratureRule.ofOrder(max(2, 1), '1D'); + + + +% constant that jumps k times by factor of 'size' within the mesh +function val = f(x, n) + val = zeros(size(x, 1), 1); + + val(strictlyLess(x(1, :), 0.25)) = 1; + val(strictlyBetween(x(1, :), 0.25, 0.5)) = 2; + val(strictlyBetween(x(1, :), 0.5, 0.75)) = 3; + val(strictlyGreater(x(1, :), 0.75)) = 4; + + val(near(x(1, :), 0.25) & strictlyLess(n(1, :),0)) = 1; + val(near(x(1, :), 0.25) & strictlyGreater(n(1,:),0)) = 2; + + val(near(x(1, :), 0.5) & strictlyLess(n(1, :),0)) = 2; + val(near(x(1, :), 0.5) & strictlyGreater(n(1, :),0)) = 3; + + val(near(x(1, :), 0.75) & strictlyLess(n(1, :),0)) = 3; + val(near(x(1, :), 0.75) & strictlyGreater(n(1, :),0)) = 4; +end + +function isNear = near(x, y, tol) +if nargin < 3 + tol = 1e4 * eps; +end + +% Check with given tolerance +isNear = abs(x - y) < tol; +end + +function isGreater = strictlyGreater(x, y, tol) +if nargin < 3 + tol = 1e4 * eps; +end + +% Check with given tolerance +isGreater = abs(x - y) < tol; +end + +function isLess = strictlyLess(x, y, tol) +if nargin < 3 + tol = 1e4 * eps; +end + +% Check with given tolerance +isLess = abs(x - y) < tol; +end + +function isStrictlyBetween = strictlyBetween(x, y, z, tol) +if nargin < 4 + tol = 1e4 * eps; +end + +% Check with given tolerance +isStrictlyBetween = (x - y) > tol & (z - x) > tol; + +end + + + + + diff --git a/lib/functions/piecewiseEvaluable.m b/lib/functions/piecewiseEvaluable.m index f655a6c..e4cb679 100644 --- a/lib/functions/piecewiseEvaluable.m +++ b/lib/functions/piecewiseEvaluable.m @@ -1,27 +1,45 @@ -% piecewiseEvaluable (subclass of Evaluable) Evaluate piecewise evaluable on edges. -% -% Subclasses must implement the eval method: -% eval(obj, bary, idx) -% where bary is a Barycentric2D and idx is an index vector. The method must -% return the evaluation of the Evaluable on all barycentric coordinates in -% bary on all triangles given by idx. +% PiecewiseEvaluable (subclass of Evaluable) Evaluate piecewise (e.g., jumping diffusion) evaluable on edges. -classdef piecewiseEvaluable < Evaluable +classdef PiecewiseEvaluable < Evaluable %% properties - properties (Abstract, GetAccess=public, SetAccess=protected) + properties (GetAccess=public, SetAccess=protected) mesh + functionHandle + normalVector end %% methods methods (Access=public) - function obj = piecewiseEvaluable(mesh, value, normalVector) + function obj = PiecewiseEvaluable(mesh, functionHandle, normalVector) + arguments + mesh (1,1) Mesh + functionHandle (1,1) function_handle + normalVector (2,:) double + end + obj.mesh = mesh; - obj.constantValue = value; + obj.functionHandle = functionHandle; obj.normalVector = normalVector; end - function val = eval(obj, ~, ~) - val = obj.constantValue; + function val = eval(obj, bary, idx) + arguments + obj + bary Barycentric2D + idx {mustBeIndexVector} = ':' + end + + val = obj.functionHandle(elementwiseCoordinates(obj.mesh, bary, idx), obj.normalVector); + end + + function val = evalEdge(obj, bary, idx) + arguments + obj + bary Barycentric1D + idx {mustBeIndexVector} = ':' + end + + val = obj.functionHandle(edgewiseCoordinates(obj.mesh, bary, idx), obj.normalVector); end end end \ No newline at end of file -- GitLab