changeset 1281:01a0500de446 feature/poroelastic

Add working implementation of Elastic2dStaggeredCurvilinearAnisotropic, which wraps around the Cartesian anisotropic scheme.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 07 Jun 2020 11:33:44 -0700
parents c7f5b480e455
children a52033540dd9 a8e730db76e9
files +scheme/Elastic2dStaggeredCurvilinearAnisotropic.m
diffstat 1 files changed, 694 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/Elastic2dStaggeredCurvilinearAnisotropic.m	Sun Jun 07 11:33:44 2020 -0700
@@ -0,0 +1,694 @@
+classdef Elastic2dStaggeredCurvilinearAnisotropic < scheme.Scheme
+
+% Discretizes the elastic wave equation:
+% rho u_{i,tt} = dj C_{ijkl} dk u_j
+% in curvilinear coordinates, using Lebedev staggered grids
+
+    properties
+        m % Number of points in each direction, possibly a vector
+        h % Grid spacing
+
+        grid
+        dim
+        nGrids
+
+        order % Order of accuracy for the approximation
+
+        % Diagonal matrices for variable coefficients
+        J, Ji
+        RHO % Density
+        C   % Elastic stiffness tensor
+
+        D  % Total operator
+
+        % Dx, Dy % Physical derivatives
+        n_w, n_e, n_s, n_n % Physical normals
+
+        % Boundary operators in cell format, used for BC
+        T_w, T_e, T_s, T_n
+
+        % Traction operators
+        % tau_w, tau_e, tau_s, tau_n      % Return vector field
+        % tau1_w, tau1_e, tau1_s, tau1_n  % Return scalar field
+        % tau2_w, tau2_e, tau2_s, tau2_n  % Return scalar field
+
+        % Inner products
+        H
+
+        % Boundary inner products (for scalar field)
+        H_w, H_e, H_s, H_n
+
+        % Surface Jacobian vectors
+        s_w, s_e, s_s, s_n
+
+        % Boundary restriction operators
+        e_w_u, e_e_u, e_s_u, e_n_u      % Act on scalar field, return scalar field at boundary
+        e_w_s, e_e_s, e_s_s, e_n_s      % Act on scalar field, return scalar field at boundary
+        % e1_w, e1_e, e1_s, e1_n  % Act on vector field, return scalar field at boundary
+        % e2_w, e2_e, e2_s, e2_n  % Act on vector field, return scalar field at boundary
+        % e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
+        % en_w, en_e, en_s, en_n  % Act on vector field, return normal component
+
+        % U{i}^T picks out component i
+        U
+
+        % G{i}^T picks out displacement grid i
+        G
+
+        % Elastic2dVariableAnisotropic object for reference domain
+        refObj
+    end
+
+    methods
+
+        % The coefficients can either be function handles or grid functions
+        % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
+        function obj = Elastic2dStaggeredCurvilinearAnisotropic(g, order, rho, C)
+            default_arg('rho', @(x,y) 0*x+1);
+
+            opSet = @sbp.D1StaggeredUpwind;
+            dim = 2;
+            nGrids = 2;
+
+            C_default = cell(dim,dim,dim,dim);
+            for i = 1:dim
+                for j = 1:dim
+                    for k = 1:dim
+                        for l = 1:dim
+                            C_default{i,j,k,l} = @(x,y) 0*x;
+                        end
+                    end
+                end
+            end
+            default_arg('C', C_default);
+            assert(isa(g, 'grid.Staggered'));
+
+            g_u = g.gridGroups{1};
+            g_s = g.gridGroups{2};
+
+            m_u = {g_u{1}.size(), g_u{2}.size()};
+            m_s = {g_s{1}.size(), g_s{2}.size()};
+
+            if isa(rho, 'function_handle')
+                rho_vec = cell(nGrids, 1);
+                for a = 1:nGrids
+                    rho_vec{a} = grid.evalOn(g_u{a}, rho);
+                end
+                rho = rho_vec;
+            end
+            for a = 1:nGrids
+                RHO{a} = spdiag(rho{a});
+            end
+            obj.RHO = RHO;
+
+            C_mat = cell(nGrids, 1);
+            for a = 1:nGrids
+                C_mat{a} = cell(dim,dim,dim,dim);
+            end
+            for a = 1:nGrids
+                for i = 1:dim
+                    for j = 1:dim
+                        for k = 1:dim
+                            for l = 1:dim
+                                if numel(C) == dim
+                                    C_mat{a}{i,j,k,l} = spdiag(C{a}{i,j,k,l});
+                                else
+                                    C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l}));
+                                end
+                            end
+                        end
+                    end
+                end
+            end
+            obj.C = C_mat;
+
+            C = cell(nGrids, 1);
+            for a = 1:nGrids
+                C{a} = cell(dim,dim,dim,dim);
+                for i = 1:dim
+                    for j = 1:dim
+                        for k = 1:dim
+                            for l = 1:dim
+                                C{a}{i,j,k,l} = diag(C_mat{a}{i,j,k,l});
+                            end
+                        end
+                    end
+                end
+            end
+
+            % Reference m for primal grid
+            m = g_u{1}.size();
+
+            % 1D operators
+            ops = cell(dim,1);
+            D1p = cell(dim, 1);
+            D1d = cell(dim, 1);
+            mp = cell(dim, 1);
+            md = cell(dim, 1);
+            Ip = cell(dim, 1);
+            Id = cell(dim, 1);
+            Hp = cell(dim, 1);
+            Hd = cell(dim, 1);
+
+            opSet = @sbp.D1StaggeredUpwind;
+            for i = 1:dim
+                ops{i} = opSet(m(i), {0,1}, order);
+                D1p{i} = ops{i}.D1_dual;
+                D1d{i} = ops{i}.D1_primal;
+                mp{i} = length(ops{i}.x_primal);
+                md{i} = length(ops{i}.x_dual);
+                Ip{i} = speye(mp{i}, mp{i});
+                Id{i} = speye(md{i}, md{i});
+                Hp{i} = ops{i}.H_primal;
+                Hd{i} = ops{i}.H_dual;
+                ep_l{i} = ops{i}.e_primal_l;
+                ep_r{i} = ops{i}.e_primal_r;
+                ed_l{i} = ops{i}.e_dual_l;
+                ed_r{i} = ops{i}.e_dual_r;
+            end
+
+            % D1_u2s{a, b}{i} approximates ddi and
+            % takes from u grid number b to s grid number a
+            % Some of D1_x2y{a, b} are 0.
+            D1_u2s = cell(nGrids, nGrids);
+            D1_s2u = cell(nGrids, nGrids);
+
+            N_u = cell(nGrids, 1);
+            N_s = cell(nGrids, 1);
+            for a = 1:nGrids
+                N_u{a} = g_u{a}.N();
+                N_s{a} = g_s{a}.N();
+            end
+
+            %---- Grid layout -------
+            % gu1 = xp o yp;
+            % gu2 = xd o yd;
+            % gs1 = xd o yp;
+            % gs2 = xp o yd;
+            %------------------------
+
+            % Logical operators
+            D1_u2s{1,1}{1} = kron(D1p{1}, Ip{2});
+            D1_s2u{1,1}{1} = kron(D1d{1}, Ip{2});
+
+            D1_u2s{1,2}{2} = kron(Id{1}, D1d{2});
+            D1_u2s{2,1}{2} = kron(Ip{1}, D1p{2});
+
+            D1_s2u{1,2}{2} = kron(Ip{1}, D1d{2});
+            D1_s2u{2,1}{2} = kron(Id{1}, D1p{2});
+
+            D1_u2s{2,2}{1} = kron(D1d{1}, Id{2});
+            D1_s2u{2,2}{1} = kron(D1p{1}, Id{2});
+
+            D1_u2s{1,1}{2} = sparse(N_s{1}, N_u{1});
+            D1_s2u{1,1}{2} = sparse(N_u{1}, N_s{1});
+
+            D1_u2s{2,2}{2} = sparse(N_s{2}, N_u{2});
+            D1_s2u{2,2}{2} = sparse(N_u{2}, N_s{2});
+
+            D1_u2s{1,2}{1} = sparse(N_s{1}, N_u{2});
+            D1_s2u{1,2}{1} = sparse(N_u{1}, N_s{2});
+
+            D1_u2s{2,1}{1} = sparse(N_s{2}, N_u{1});
+            D1_s2u{2,1}{1} = sparse(N_u{2}, N_s{1});
+
+
+            %---- Combine grids and components -----
+
+            % U{a}{i}^T picks out u component i on grid a
+            U = cell(nGrids, 1);
+            for a = 1:nGrids
+                U{a} = cell(dim, 1);
+                I = speye(N_u{a}, N_u{a});
+                for i = 1:dim
+                    E = sparse(dim,1);
+                    E(i) = 1;
+                    U{a}{i} = kron(I, E);
+                end
+            end
+            obj.U = U;
+
+            % Order grids
+            % u1, u2
+            Iu1 = speye(dim*N_u{1}, dim*N_u{1});
+            Iu2 = speye(dim*N_u{2}, dim*N_u{2});
+
+            Gu1 = cell2mat( {Iu1; sparse(dim*N_u{2}, dim*N_u{1})} );
+            Gu2 = cell2mat( {sparse(dim*N_u{1}, dim*N_u{2}); Iu2} );
+
+            G = {Gu1; Gu2};
+
+            %---- Grid layout -------
+            % gu1 = xp o yp;
+            % gu2 = xd o yd;
+            % gs1 = xd o yp;
+            % gs2 = xp o yd;
+            %------------------------
+
+            % Boundary restriction ops
+            e_w_u = cell(nGrids, 1);
+            e_s_u = cell(nGrids, 1);
+            e_e_u = cell(nGrids, 1);
+            e_n_u = cell(nGrids, 1);
+
+            e_w_s = cell(nGrids, 1);
+            e_s_s = cell(nGrids, 1);
+            e_e_s = cell(nGrids, 1);
+            e_n_s = cell(nGrids, 1);
+
+            e_w_u{1} = kron(ep_l{1}, Ip{2});
+            e_e_u{1} = kron(ep_r{1}, Ip{2});
+            e_s_u{1} = kron(Ip{1}, ep_l{2});
+            e_n_u{1} = kron(Ip{1}, ep_r{2});
+
+            e_w_u{2} = kron(ed_l{1}, Id{2});
+            e_e_u{2} = kron(ed_r{1}, Id{2});
+            e_s_u{2} = kron(Id{1}, ed_l{2});
+            e_n_u{2} = kron(Id{1}, ed_r{2});
+
+            e_w_s{1} = kron(ed_l{1}, Ip{2});
+            e_e_s{1} = kron(ed_r{1}, Ip{2});
+            e_s_s{1} = kron(Id{1}, ep_l{2});
+            e_n_s{1} = kron(Id{1}, ep_r{2});
+
+            e_w_s{2} = kron(ep_l{1}, Id{2});
+            e_e_s{2} = kron(ep_r{1}, Id{2});
+            e_s_s{2} = kron(Ip{1}, ed_l{2});
+            e_n_s{2} = kron(Ip{1}, ed_r{2});
+
+            %  --- Metric coefficients on stress grids -------
+            x = cell(nGrids, 1);
+            y = cell(nGrids, 1);
+            J = cell(nGrids, 1);
+            x_xi = cell(nGrids, 1);
+            x_eta = cell(nGrids, 1);
+            y_xi = cell(nGrids, 1);
+            y_eta = cell(nGrids, 1);
+
+            for a = 1:nGrids
+                coords = g_u{a}.points();
+                x{a} = coords(:,1);
+                y{a} = coords(:,2);
+            end
+
+            for a = 1:nGrids
+                x_xi{a} = zeros(N_s{a}, 1);
+                y_xi{a} = zeros(N_s{a}, 1);
+                x_eta{a} = zeros(N_s{a}, 1);
+                y_eta{a} = zeros(N_s{a}, 1);
+
+                for b = 1:nGrids
+                    x_xi{a} = x_xi{a} + D1_u2s{a,b}{1}*x{b};
+                    y_xi{a} = y_xi{a} + D1_u2s{a,b}{1}*y{b};
+                    x_eta{a} = x_eta{a} + D1_u2s{a,b}{2}*x{b};
+                    y_eta{a} = y_eta{a} + D1_u2s{a,b}{2}*y{b};
+                end
+            end
+
+            for a = 1:nGrids
+                J{a} = x_xi{a}.*y_eta{a} - x_eta{a}.*y_xi{a};
+            end
+
+            K = cell(nGrids, 1);
+            for a = 1:nGrids
+                K{a} = cell(dim, dim);
+                K{a}{1,1} = y_eta{a}./J{a};
+                K{a}{1,2} = -y_xi{a}./J{a};
+                K{a}{2,1} = -x_eta{a}./J{a};
+                K{a}{2,2} = x_xi{a}./J{a};
+            end
+            % ----------------------------------------------
+
+            %  --- Metric coefficients on displacement grids -------
+            x_s = cell(nGrids, 1);
+            y_s = cell(nGrids, 1);
+            J_u = cell(nGrids, 1);
+            x_xi_u = cell(nGrids, 1);
+            x_eta_u = cell(nGrids, 1);
+            y_xi_u = cell(nGrids, 1);
+            y_eta_u = cell(nGrids, 1);
+
+            for a = 1:nGrids
+                coords = g_s{a}.points();
+                x_s{a} = coords(:,1);
+                y_s{a} = coords(:,2);
+            end
+
+            for a = 1:nGrids
+                x_xi_u{a} = zeros(N_u{a}, 1);
+                y_xi_u{a} = zeros(N_u{a}, 1);
+                x_eta_u{a} = zeros(N_u{a}, 1);
+                y_eta_u{a} = zeros(N_u{a}, 1);
+
+                for b = 1:nGrids
+                    x_xi_u{a} = x_xi_u{a} + D1_s2u{a,b}{1}*x_s{b};
+                    y_xi_u{a} = y_xi_u{a} + D1_s2u{a,b}{1}*y_s{b};
+                    x_eta_u{a} = x_eta_u{a} + D1_s2u{a,b}{2}*x_s{b};
+                    y_eta_u{a} = y_eta_u{a} + D1_s2u{a,b}{2}*y_s{b};
+                end
+            end
+
+            for a = 1:nGrids
+                J_u{a} = x_xi_u{a}.*y_eta_u{a} - x_eta_u{a}.*y_xi_u{a};
+            end
+            % ----------------------------------------------
+
+            % x_u = Du*x;
+            % x_v = Dv*x;
+            % y_u = Du*y;
+            % y_v = Dv*y;
+
+            % J = x_u.*y_v - x_v.*y_u;
+
+            % K = cell(dim, dim);
+            % K{1,1} = y_v./J;
+            % K{1,2} = -y_u./J;
+            % K{2,1} = -x_v./J;
+            % K{2,2} = x_u./J;
+
+            % Physical derivatives
+            % obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
+            % obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
+
+            % Wrap around Aniosotropic Cartesian. Transformed density and stiffness
+            rho_tilde = cell(nGrids, 1);
+            PHI = cell(nGrids, 1);
+
+            for a = 1:nGrids
+                rho_tilde{a} = J_u{a}.*rho{a};
+            end
+
+            for a = 1:nGrids
+                PHI{a} = cell(dim,dim,dim,dim);
+                for i = 1:dim
+                    for j = 1:dim
+                        for k = 1:dim
+                            for l = 1:dim
+                                PHI{a}{i,j,k,l} = 0*C{a}{i,j,k,l};
+                                for m = 1:dim
+                                    for n = 1:dim
+                                        PHI{a}{i,j,k,l} = PHI{a}{i,j,k,l} + J{a}.*K{a}{m,i}.*C{a}{m,j,n,l}.*K{a}{n,k};
+                                    end
+                                end
+                            end
+                        end
+                    end
+                end
+            end
+
+            refObj = scheme.Elastic2dStaggeredAnisotropic(g.logic, order, rho_tilde, PHI);
+
+            G = refObj.G;
+            U = refObj.U;
+            H_u = refObj.H_u;
+
+            % Volume quadrature
+            [m, n] = size(refObj.H);
+            obj.H = sparse(m, n);
+            obj.J = sparse(m, n);
+            for a = 1:nGrids
+                for i = 1:dim
+                    obj.H = obj.H + G{a}*U{a}{i}*spdiag(J_u{a})*refObj.H_u{a}*U{a}{i}'*G{a}';
+                    obj.J = obj.J + G{a}*U{a}{i}*spdiag(J_u{a})*U{a}{i}'*G{a}';
+                end
+            end
+            obj.Ji = inv(obj.J);
+
+            % Boundary quadratures on stress grids
+            s_w = cell(nGrids, 1);
+            s_e = cell(nGrids, 1);
+            s_s = cell(nGrids, 1);
+            s_n = cell(nGrids, 1);
+
+            % e_w_u = refObj.e_w_u;
+            % e_e_u = refObj.e_e_u;
+            % e_s_u = refObj.e_s_u;
+            % e_n_u = refObj.e_n_u;
+
+            e_w_s = refObj.e_w_s;
+            e_e_s = refObj.e_e_s;
+            e_s_s = refObj.e_s_s;
+            e_n_s = refObj.e_n_s;
+
+            for a = 1:nGrids
+                s_w{a} = sqrt((e_w_s{a}'*x_eta{a}).^2 + (e_w_s{a}'*y_eta{a}).^2);
+                s_e{a} = sqrt((e_e_s{a}'*x_eta{a}).^2 + (e_e_s{a}'*y_eta{a}).^2);
+                s_s{a} = sqrt((e_s_s{a}'*x_xi{a}).^2 + (e_s_s{a}'*y_xi{a}).^2);
+                s_n{a} = sqrt((e_n_s{a}'*x_xi{a}).^2 + (e_n_s{a}'*y_xi{a}).^2);
+            end
+
+            obj.s_w = s_w;
+            obj.s_e = s_e;
+            obj.s_s = s_s;
+            obj.s_n = s_n;
+
+            % for a = 1:nGrids
+                % obj.H_w_s{a} = refObj.H_w_s{a}*spdiag(s_w{a});
+                % obj.H_e_s{a} = refObj.H_e_s{a}*spdiag(s_e{a});
+                % obj.H_s_s{a} = refObj.H_s_s{a}*spdiag(s_s{a});
+                % obj.H_n_s{a} = refObj.H_n_s{a}*spdiag(s_n{a});
+            % end
+
+            % Restriction operators
+            obj.e_w_u = refObj.e_w_u;
+            obj.e_e_u = refObj.e_e_u;
+            obj.e_s_u = refObj.e_s_u;
+            obj.e_n_u = refObj.e_n_u;
+
+            obj.e_w_s = refObj.e_w_s;
+            obj.e_e_s = refObj.e_e_s;
+            obj.e_s_s = refObj.e_s_s;
+            obj.e_n_s = refObj.e_n_s;
+
+            % Adapt things from reference object
+            obj.D = refObj.D;
+            obj.U = refObj.U;
+            obj.G = refObj.G;
+
+            % obj.e1_w = refObj.e1_w;
+            % obj.e1_e = refObj.e1_e;
+            % obj.e1_s = refObj.e1_s;
+            % obj.e1_n = refObj.e1_n;
+
+            % obj.e2_w = refObj.e2_w;
+            % obj.e2_e = refObj.e2_e;
+            % obj.e2_s = refObj.e2_s;
+            % obj.e2_n = refObj.e2_n;
+
+            % obj.e_scalar_w = refObj.e_scalar_w;
+            % obj.e_scalar_e = refObj.e_scalar_e;
+            % obj.e_scalar_s = refObj.e_scalar_s;
+            % obj.e_scalar_n = refObj.e_scalar_n;
+
+            % e1_w = obj.e1_w;
+            % e1_e = obj.e1_e;
+            % e1_s = obj.e1_s;
+            % e1_n = obj.e1_n;
+
+            % e2_w = obj.e2_w;
+            % e2_e = obj.e2_e;
+            % e2_s = obj.e2_s;
+            % e2_n = obj.e2_n;
+
+            % obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
+            % obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
+            % obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
+            % obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
+
+            % obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
+            % obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
+            % obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
+            % obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
+
+            % obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
+            % obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
+            % obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
+            % obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
+
+            % % Physical normals
+            % e_w = obj.e_scalar_w;
+            % e_e = obj.e_scalar_e;
+            % e_s = obj.e_scalar_s;
+            % e_n = obj.e_scalar_n;
+
+            % e_w_vec = obj.e_w;
+            % e_e_vec = obj.e_e;
+            % e_s_vec = obj.e_s;
+            % e_n_vec = obj.e_n;
+
+            % nu_w = [-1,0];
+            % nu_e = [1,0];
+            % nu_s = [0,-1];
+            % nu_n = [0,1];
+
+            % obj.n_w = cell(2,1);
+            % obj.n_e = cell(2,1);
+            % obj.n_s = cell(2,1);
+            % obj.n_n = cell(2,1);
+
+            % n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
+            % n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
+            % obj.n_w{1} = spdiag(n_w_1);
+            % obj.n_w{2} = spdiag(n_w_2);
+
+            % n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
+            % n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
+            % obj.n_e{1} = spdiag(n_e_1);
+            % obj.n_e{2} = spdiag(n_e_2);
+
+            % n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
+            % n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
+            % obj.n_s{1} = spdiag(n_s_1);
+            % obj.n_s{2} = spdiag(n_s_2);
+
+            % n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
+            % n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
+            % obj.n_n{1} = spdiag(n_n_1);
+            % obj.n_n{2} = spdiag(n_n_2);
+
+            % % Operators that extract the normal component
+            % obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
+            % obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
+            % obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
+            % obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
+
+            % Misc.
+            obj.refObj = refObj;
+            obj.m = refObj.m;
+            obj.h = refObj.h;
+            obj.order = order;
+            obj.grid = g;
+            obj.dim = dim;
+            obj.nGrids = nGrids;
+
+        end
+
+
+        % Closure functions return the operators applied to the own domain to close the boundary
+        % Penalty functions return the operators 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'.
+        %       bc                  is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
+        %                           on the first component. Can also be e.g.
+        %                           {'normal', 'd'} or {'tangential', 't'} for conditions on
+        %                           tangential/normal component.
+        %       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.
+
+        % For displacement bc:
+        % bc = {comp, 'd', dComps},
+        % where
+        % dComps = vector of components with displacement BC. Default: 1:dim.
+        % In this way, we can specify one BC at a time even though the SATs depend on all BC.
+        function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
+            default_arg('tuning', 1.0);
+            assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
+
+            [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
+
+            type = bc{2};
+
+            switch type
+            case {'F','f','Free','free','traction','Traction','t','T'}
+                s = obj.(['s_' boundary]);
+                s = spdiag(cell2mat(s));
+                penalty = penalty*s;
+            end
+        end
+
+        % type     Struct that specifies the interface coupling.
+        %          Fields:
+        %          -- tuning:           penalty strength, defaults to 1.0
+        %          -- interpolation:    type of interpolation, default 'none'
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+
+            defaultType.tuning = 1.0;
+            defaultType.interpolation = 'none';
+            default_struct('type', defaultType);
+
+            [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
+        end
+
+        % Returns h11 for the boundary specified by the string boundary.
+        % op -- string
+        function h11 = getBorrowing(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
+            switch boundary
+            case {'w','e'}
+                h11 = obj.refObj.h11{1};
+            case {'s', 'n'}
+                h11 = obj.refObj.h11{2};
+            end
+        end
+
+        % Returns the outward unit normal vector for the boundary specified by the string boundary.
+        % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
+        function n = getNormal(obj, boundary)
+            error('Not implemented');
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
+            n = obj.(['n_' boundary]);
+        end
+
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op -- string
+        function o = getBoundaryOperator(obj, op, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
+
+            o = obj.([op, '_', boundary]);
+
+        end
+
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op -- string
+        function o = getBoundaryOperatorForScalarField(obj, op, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+            assertIsMember(op, {'e_u', 'e_s'})
+
+            switch op
+                case 'e_u'
+                    o = obj.(['e_', boundary, '_u']);
+                case 'e_s'
+                    o = obj.(['e_', boundary, '_s']);
+            end
+
+        end
+
+        % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
+        % Formula: tau_i = T_ij u_j
+        % op -- string
+        function T = getBoundaryTractionOperator(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
+            T = obj.(['T', '_', boundary]);
+        end
+
+        % Returns square boundary quadrature matrix, of dimension
+        % corresponding to the number of boundary unknowns
+        %
+        % boundary -- string
+        function H = getBoundaryQuadrature(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
+            H = obj.getBoundaryQuadratureForScalarField(boundary);
+            I_dim = speye(obj.dim, obj.dim);
+            H = kron(H, I_dim);
+        end
+
+        % Returns square boundary quadrature matrix, of dimension
+        % corresponding to the number of boundary grid points
+        %
+        % boundary -- string
+        function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
+            H_b = obj.(['H_', boundary]);
+        end
+
+        function N = size(obj)
+            N = length(obj.D);
+        end
+    end
+end