Mercurial > repos > public > sbplib
view +scheme/Laplace1dVariable.m @ 1324:56439c1a49b4 feature/laplace1d_variable
Bugifx penalty sign.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Sun, 27 Jun 2021 16:19:28 +0200 |
parents | 3f4e011193f0 |
children |
line wrap: on
line source
classdef Laplace1dVariable < scheme.Scheme properties grid order % Order accuracy for the approximation D % scheme operator H % Discrete norm % Variable coefficients a(b u_x)_x a b Hi e_l e_r d_l d_r gamm end methods function obj = Laplace1dVariable(g, order, a, b) default_arg('a', @(x) 0*x + 1); default_arg('b', @(x) 0*x + 1); assertType(g, 'grid.Cartesian'); ops = sbp.D2Variable(g.size(), g.lim{1}, order); obj.H = sparse(ops.H); obj.Hi = sparse(ops.HI); obj.e_l = sparse(ops.e_l); obj.e_r = sparse(ops.e_r); obj.d_l = -sparse(ops.d1_l); obj.d_r = sparse(ops.d1_r); obj.grid = g; obj.order = order; a = grid.evalOn(g, a); b = grid.evalOn(g, b); A = spdiag(a); B = spdiag(b); obj.a = A; obj.b = B; obj.D = A*ops.D2(b); obj.gamm = g.h*ops.borrowing.M.d1; end % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. % type is a string specifying the type of boundary condition if there are several. % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj,boundary,type,data) default_arg('type','neumann'); default_arg('data',0); e = obj.getBoundaryOperator('e', boundary); d = obj.getBoundaryOperator('d', boundary); s = obj.getBoundarySign(boundary); b = e'*obj.b*e; switch type % Dirichlet boundary condition case {'D','d','dirichlet'} tuning = 1.1; tau1 = -tuning/obj.gamm; tau2 = 1; tau = tau1*e + tau2*d; closure = obj.a*obj.Hi*tau*b*e'; penalty = -obj.a*obj.Hi*tau*b; % Neumann boundary condition case {'N','n','neumann'} tau = -e; closure = obj.a*obj.Hi*tau*b*d'; penalty = -obj.a*obj.Hi*tau*b; % Unknown, boundary condition otherwise error('No such boundary condition: type = %s',type); end end function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain e_u = obj.getBoundaryOperator('e', boundary); d_u = obj.getBoundaryOperator('d', boundary); s_u = obj.getBoundarySign(boundary); e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary); s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); b_u = e_u'*obj.b*e_u; b_v = e_v'*neighbour_scheme.b*e_v; gamm_u = obj.gamm; gamm_v = neighbour_scheme.gamm; tuning = 1.1; closure = zeros(size(obj.D)); penalty = zeros(length(obj.D), length(b_v)); % Continuity of bu_x closure = closure - 1/2*e_u*b_u*d_u'; penalty = penalty - 1/2*e_u*b_v*d_v'; % Continuity of u (symmetrizing term) closure = closure + 1/2*d_u*b_u*e_u'; penalty = penalty - 1/2*d_u*b_u*e_v'; % Continuity of u (symmetric term) tau = 1/4*(b_u/gamm_u + b_v/gamm_v)*tuning; closure = closure - e_u*tau*e_u'; penalty = penalty + e_u*tau*e_v'; % Multiply by Hi and a. closure = obj.Hi*obj.a*closure; penalty = obj.Hi*obj.a*penalty; end % Returns the boundary operator op for the boundary specified by the string boundary. % op -- string % boundary -- string function o = getBoundaryOperator(obj, op, boundary) assertIsMember(op, {'e', 'd'}) assertIsMember(boundary, {'l', 'r'}) o = obj.([op, '_', boundary]); end % Returns square boundary quadrature matrix, of dimension % corresponding to the number of boundary points % % boundary -- string % Note: for 1d diffOps, the boundary quadrature is the scalar 1. function H_b = getBoundaryQuadrature(obj, boundary) assertIsMember(boundary, {'l', 'r'}) H_b = 1; end % Returns the boundary sign. The right boundary is considered the positive boundary % boundary -- string function s = getBoundarySign(obj, boundary) assertIsMember(boundary, {'l', 'r'}) switch boundary case {'r'} s = 1; case {'l'} s = -1; end end function N = size(obj) N = obj.grid.size(); end end end