FEM solver for electrostatics report. Sauter and schwab quadrature and required tools

Forcecomputation.m

+% Size of the mesh
+N = 100;
+
+% Dimensions of mesh?
+L = [1 1 1];
+
+% Creating a mesh
+m = mshCube(N,L);
+
+% Defining integration domains?
+omega = dom(m,4); % 4 gauss points
+
+% Example of manipulating a mesh
+m2 = mshCube(N,2*L);
+m2 = translate(m2, [5 5 5]);
+
+% Plotting a mest
+plot(m2);
+
+% Getting the boundary of the mesh as a mesh object
+mbound = mshBoundary(m);
+
+% domain of integration
+Gamma = dom(mbound,3);
+
+% Creating the spaces
+
+% Piecewise linear functions
+Vh = fem(m,'P1');
+
+% Ravier thomas 'RWG', Nedelec 'Ned' for other types 
+
+
+% A way to specify the Dirichlet conditions? Vh had an element called dir?
+Vh0 = dirichlet(Vh, mbound);
+
+% Matlab anonymous function? X has 3 columns, function returns the values
+% vectorially
+G = @(X) (sin(X(:,1)) .* cos(X(:,2)) .* X(:,3));
+
+% Creating a function for grad G?
+nablaG = cell(3,1);
+nablaG{1} = @(X)(cos(X(:,1).*cos(X(:,2)) .*)s
+
+% Creating integrals 
+rhs = integral(Omega,grad(Vh0),gradG);
+
+% Bilinear form
+a = integral(Omega, grad(Vh0),grad(Vh0));
+
+%solution
+vh = a\rhs;
+
+% Creating a finite element restriction for G
+%Gh = G(m.vtx); % does not work because we have to remove boundary points?
+
+% m.vtx contains the boundary that's why we do it this way
+%Gh = G(Vh0.dof); % dof for p2 will also contain centers
+Gh = G(Vh0.unk); % What is this?
+
+% What is the feval function?
+
+
+

Gypsi_files/addpathGypsilab.m

+function addpathGypsilab()
+addpath('miscellaneous')
+addpath('openBmm')
+addpath('openDom')
+addpath('openEbd')
+addpath('openFem')
+addpath('openFfm')
+addpath('openHmx')
+addpath('openMmg')
+addpath('openMsh')
+addpath('openOpr')
+addpath('openRay')
+end
\ No newline at end of file

Gypsi_files/openFem/fem.m

@@ -255,6 +255,11 @@ methods
         M = femUnk2Qud(fe,domain);
     end
     
+    % Get the reference shape functions
+    function rsfs = rsf(fe)
+       rsfs = femRSF(fe); 
+    end
+    
     
     
     % UNKNOWNS DATA TO VERTEX

Gypsi_files/openFem/femRSF.m

+%% Getting the hard coded reference shape functions for the FEM object
+
+% IMPORTANT: The triangular reference element is different than the one
+% used in Gypsilab.
+
+function rsfs = femRSF(fe)
+    % Getting the type of FE space
+    typ = fe.typ;
+    mesh = fe.msh;
+    
+    elt_type = size(mesh.elt,2);
+    
+    % Segment (Is there a reference element in gypsi for line segment?)
+    if elt_type == 2
+        switch typ
+            case 'P0'
+                rsfs = cell(1,1);
+                rsfs{1} = @(X) 1;
+            case 'P1'
+                rsfs = cell(2,1);
+                rsfs{1} = @(X) 0.5 * (1+X); % support at 1
+                rsfs{2} = @(X) 0.5 * (1-X); % support at -1
+        end
+        
+    end
+    % Triangle
+    if elt_type == 3
+       switch typ
+            case 'P0'
+                rsfs = cell(1,1);
+                rsfs{1} = @(X) 1;
+            case 'P1'
+                rsfs = cell(3,1);
+                %rsfs{1} = @(X) 1 - X(1) - X(2);
+                %rsfs{2} = @(X) X(1);
+                %rsfs{3} = @(X) X(2);
+                rsfs{1} = @(X) 1 - X(1); % support at 0,0
+                rsfs{2} = @(X) X(1) - X(2); % support at 1,0 
+                rsfs{3} = @(X) X(2); % support at 1,1
+       end
+    end
+end
\ No newline at end of file

Gypsi_files/openMsh/mshReadMesh.m

@@ -37,8 +37,11 @@ end
 
 % Nodes
 str = fgets(fid);
+i=0;
 while isempty(strfind(str,'Vertices'))
     str = fgets(fid);
+    i = i+1;
+    disp(i);
 end
 Nvtx = str2double(fgets(fid));
 vtx  = zeros(Nvtx,3);

Piyush's codes/Dirichlet_BVP_BEM.m

+% Solving a Dirichlet BVP on a circle
+
+close all;
+clear all;
+
+%% Creating a boundary mesh for a 2D circular domain
+
+R = 1.2; % Radius of the circle
+% Center assumed to be the origin
+
+% Mesh parameters
+N = 50;
+
+% Creating the mesh manually by constructing the vertices, elements
+theta = linspace(0, 2*pi, N+1);
+theta = theta(2:N+1);
+
+vtx = [ R * cos(theta); R * sin(theta); 0 * theta]'; % Vertices
+elt = [1:N; [2:N] 1]'; % Elements
+col = zeros(N,1); % Col?
+
+mesh = msh(vtx,elt,col);
+plot(mesh);
+
+% This can also be done by making a 2D mesh and then extracting the
+% boundary out of it
+% meshvol = mshDisk(N,R);
+% mesh = meshvol.bnd;
+
+%% Setting up the discrete system
+
+% Quadrature order
+order = 3;
+% Creating the integration domain
+Gamma = dom(mesh, order);
+
+% Specifying the boundary conditions using function g
+%g = @(X)sin(X(:,1)) .* cos(X(:,2));
+%g = @(X)ones(size(X(:,1)))*0;
+g = @(X)sin(X(:,1)-X(:,2)) .* sinh(X(:,1)+X(:,2));
+
+% Defining the FEM spaces
+trial_space = fem(mesh,'P1');
+test_space = fem(mesh,'P1');
+% A different trial space for the boundary condition g?
+trial_space_g = fem(mesh,'P1');
+
+% Defining the kernel for SL
+Gxy = @(X,Y)femGreenKernel(X,Y,'[log(r)]',0); % 0 wave number
+% Defining the kernel for DL
+GradGn = cell(3,1);
+GradGn{1} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]1',0);
+GradGn{2} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]2',0);
+GradGn{3} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]3',0);
+
+% LHS
+% Evaluating the bilinear form for the Single Layer BIO
+V = -1/(2*pi)*integral(Gamma,Gamma,trial_space,Gxy,test_space);
+V = V + -1/(2*pi)*regularize(Gamma,Gamma,trial_space,'[log(r)]',test_space);
+
+% RHS
+% Evaluating the bilinear form for the Double Layer BIO
+K = -1/(2*pi)*integral(Gamma,Gamma,trial_space,GradGn,ntimes(test_space));
+% Regularization
+K = K -1/(2*pi)*regularize(Gamma,Gamma,trial_space,'grady[log(r)]',ntimes(test_space));
+
+l = integral(Gamma,trial_space,test_space);
+
+% A clean way to get the coefficients g_N such that when combined with the
+% basis elements in trial_space_g, we get an approximation of g?
+g_N = g(trial_space_g.dof);
+
+% The linear system is V x = (0.5 l + K) g_N;
+
+% Neumann trace solution:
+Psi = V \ ((0.5 * l + K) * g_N );
+
+%% Computing the exact solution for testing accuracy
+
+x = R * cos(theta);
+y = R * sin(theta);
+
+f_x = cos(x-y) .* sinh(x+y) + sin(x-y) .* cosh(x+y);
+f_y = -cos(x-y) .* sinh(x+y) + sin(x-y) .* cosh(x+y);
+
+Psi_exact = f_x .* cos(theta) + f_y .* sin(theta);
+
+figure;
+
+plot(theta,Psi);
+hold on;
+plot(theta,Psi_exact);
+
+% Exact solution matches with computed solution!
+
+%% Computing without discretizing the boundary condition g
+% LHS remains the same, computing RHS:
+l_ex = integral(Gamma,trial_space,g);
+
+K_ex = -1/(2*pi)*integral(Gamma,Gamma,g,GradGn,ntimes(test_space));
+% Regularization
+K_ex = K_ex -1/(2*pi)*regularize(Gamma,Gamma,g,'grady[log(r)]',ntimes(test_space));
+
+% Neumann trace solution:
+Psi_ex = V \ (0.5 * l_ex + K_ex);
+
+
+%% Solving the adjoint equation
+
+% The LHS uses the same V matrix computed above. Computing the RHS:
+
+RHS_adj = 0.5 ;
+
+%% Martin stuff
+
+% uqm is unknown to quadrature matrix, used on a fem object
+% qud is the quadrature points?
+
+
+
+
+
+
+
+

Piyush's codes/Dirichlet_BVP_FEM.m

+% Solving a Dirichlet BVP on a circle
+
+close all;
+clear all;
+%opengl('save','software')
+% Recursive way by genpath? LOL
+%addpath(genpath('pathtothatfolder')) 
+% ! -> allows to read bash
+% e.g. !cd ..
+%% Creating a volume mesh for a 2D circular domain
+
+R = 1.2; % Radius of the circle centered at origin
+
+% Mesh parameters
+N = 4;
+
+% Using the in built function to create the mesh
+mesh = mshDisk(N,R);
+bnd_mesh = mesh.bnd;
+plot(mesh);
+
+%% Setting up the discrete system
+
+% Quadrature order
+order = 3;
+
+% Creating the integration domain
+Omega = dom(mesh, order);
+
+% Specifying the boundary conditions using function g
+%g = @(X)sin(X(:,1)) .* cos(X(:,2));
+%g = @(X)ones(size(X(:,1)))*0;
+g = @(X)sin(X(:,1)-X(:,2)) .* sinh(X(:,1)+X(:,2));
+
+% FEM space without the dirichlet boundary condition
+space = fem(mesh,'P1');
+% FEM space with the homogeneous Dirichlet BC
+space0 = dirichlet(space, bnd_mesh);
+
+% LHS matrix
+A = integral(Omega,grad(space0),grad(space0));
+
+% How to get the RHS vector? assuming f = -Delta G
+% v = -integral(Gamma,f,space0); % Instead of this, use integration by
+% parts to convert this expression with only first derivatives of G
+
+% Is it necessary to discretize G? using the extension by zero trick!
+% Evaluating g at mesh.vtx would require knowledge of the function on the
+% interior so we need zeros at the positions where vertices are in the 
+% interior of the mesh
+
+% Initializing to a vector of ones and zeros
+%G_N = zeros(size(space.unk,1),1);
+%
+u_ex = g(space.unk);
+
+% Naive approach
+%idx = find(abs(vecnorm(space.unk,2,2)-R) < 1e-12);
+%P = zeros(size(space.unk,1),1);
+%P(idx) = 1;
+%G_N = u_ex .* P;
+
+% Martin's restriction operator approach
+P = restriction(space,bnd_mesh); % Gives a matrix with size N_bnd X N_space
+g_N = g(bnd_mesh.vtx);
+G_N = P' * g_N; % Extension by zero
+
+v = -integral(Omega,grad(space0),grad(space))* G_N;
+
+% This is the offset solution u* 
+sol = A\v; 
+% Sol contains N_space - N_bnd elements, how to get a solution of size N_space?
+% Using Martin's approach of "elimination"
+Q = elimination(space,bnd_mesh); % Matrix of size N_space X (N_space-N_bnd)
+
+% Solution with full size N_space, extension by zero
+sol = Q * sol;
+
+% Obtaining the solution u by undoing the offset
+u = sol + G_N;
+
+err = norm(u-u_ex)
+
+%% Computing the exact solution for testing accuracy
+
+x = R * cos(theta);
+y = R * sin(theta);
+
+f_x = cos(x-y) .* sinh(x+y) + sin(x-y) .* cosh(x+y);
+f_y = -cos(x-y) .* sinh(x+y) + sin(x-y) .* cosh(x+y);
+
+Psi_exact = f_x .* cos(theta) + f_y .* sin(theta);
+
+figure;
+
+plot(theta,Psi);
+hold on;
+plot(theta,Psi_exact);
+
+% Exact solution matches with computed solution!
+
+%% Computing without discretizing the boundary condition g
+% LHS remains the same, computing RHS:
+l_ex = integral(Omega,trial_space,g);
+
+K_ex = -1/(2*pi)*integral(Omega,Omega,g,GradGn,ntimes(test_space));
+% Regularization
+K_ex = K_ex -1/(2*pi)*regularize(Omega,Omega,g,'grady[log(r)]',ntimes(test_space));
+
+% Neumann trace solution:
+Psi_ex = V \ (0.5 * l_ex + K_ex);
+
+
+%% Solving the adjoint equation
+
+% The LHS uses the same V matrix computed above. Computing the RHS:
+
+RHS_adj = 0.5 ;
+
+
+
+
+
+
+
+
+
+

Piyush's codes/Dirichlet_DFK.m

+% This function is used to solve for the Neumann Trace using the Direct 
+% First Kind BIE. Inputs: boundary mesh, function for boundary data g 
+
+function psi = Dirichlet_DFK(mesh,g)
+% Setting up the discrete system
+
+% Quadrature order
+order = 3;
+% Creating the integration domain
+Gamma = dom(mesh, order);
+
+% Defining the BEM spaces
+S1_Gamma = fem(mesh,'P1'); % Space for discretizing g
+S0_Gamma = fem(mesh,'P0'); % Usual space for Neumann trace space?
+
+% Defining the kernel for SL
+Gxy = @(X,Y)femGreenKernel(X,Y,'[log(r)]',0); % 0 wave number
+% Defining the kernel for DL
+GradGn = cell(3,1);
+GradGn{1} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]1',0);
+GradGn{2} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]2',0);
+GradGn{3} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]3',0);
+
+% LHS
+% Evaluating the bilinear form for the Single Layer BIO
+%V = -1/(2*pi)*integral(Gamma,Gamma,S1_Gamma,Gxy,S1_Gamma);
+%V = V + -1/(2*pi)*regularize(Gamma,Gamma,S1_Gamma,'[log(r)]',S1_Gamma);
+
+V = -1/(2*pi)*integral(Gamma,Gamma,S0_Gamma,Gxy,S0_Gamma);
+V = V + -1/(2*pi)*regularize(Gamma,Gamma,S0_Gamma,'[log(r)]',S0_Gamma);
+
+% RHS
+% Evaluating the bilinear form for the Double Layer BIO
+K = -1/(2*pi)*integral(Gamma,Gamma,S0_Gamma,GradGn,ntimes(S1_Gamma));
+% Regularization
+K = K -1/(2*pi)*regularize(Gamma,Gamma,S0_Gamma,'grady[log(r)]',ntimes(S1_Gamma));
+
+l = integral(Gamma,S0_Gamma,S1_Gamma);
+
+% A clean way to get the coefficients g_N such that when combined with the
+% basis elements in trial_space_g, we get an approximation of g?
+g_N = g(S1_Gamma.dof);
+
+% The linear system is V x = (0.5 l + K) g_N;
+
+% Neumann trace solution:
+psi = V \ ((0.5 * l + K) * g_N );
+    
+end
\ No newline at end of file

Piyush's codes/Dirichlet_DFK.m~

+% This function is used to solve for the Neumann Trace using the Direct 
+% First Kind BIE. Inputs: boundary mesh, function for boundary data g, 
+
+function psi = solve_dirichlet_dfk(bnd_mesh,g)
+%% Setting up the discrete system
+
+% Quadrature order
+order = 3;
+% Creating the integration domain
+Gamma = dom(mesh, order);
+
+% Defining the FEM spaces
+trial_space = fem(mesh,'P1');
+test_space = fem(mesh,'P1');
+% A different trial space for the boundary condition g?
+trial_space_g = fem(mesh,'P1');
+
+% Defining the kernel for SL
+Gxy = @(X,Y)femGreenKernel(X,Y,'[log(r)]',0); % 0 wave number
+% Defining the kernel for DL
+GradGn = cell(3,1);
+GradGn{1} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]1',0);
+GradGn{2} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]2',0);
+GradGn{3} = @(X,Y)femGreenKernel(X,Y,'grady[log(r)]3',0);
+
+% LHS
+% Evaluating the bilinear form for the Single Layer BIO
+V = -1/(2*pi)*integral(Gamma,Gamma,trial_space,Gxy,test_space);
+V = V + -1/(2*pi)*regularize(Gamma,Gamma,trial_space,'[log(r)]',test_space);
+
+% RHS
+% Evaluating the bilinear form for the Double Layer BIO
+K = -1/(2*pi)*integral(Gamma,Gamma,trial_space,GradGn,ntimes(test_space));
+% Regularization
+K = K -1/(2*pi)*regularize(Gamma,Gamma,trial_space,'grady[log(r)]',ntimes(test_space));
+
+l = integral(Gamma,trial_space,test_space);
+
+% A clean way to get the coefficients g_N such that when combined with the
+% basis elements in trial_space_g, we get an approximation of g?
+g_N = g(trial_space_g.dof);
+
+% The linear system is V x = (0.5 l + K) g_N;
+
+% Neumann trace solution:
+Psi = V \ ((0.5 * l + K) * g_N );
+    
+end
+
+
+
+clear all;
+close all;
+
+%% msh Mesh 
+mesh = msh('sqsq8.msh');
+figure;
+title('Volume mesh');
+plot(mesh);
+figure;
+bnd_mesh = mesh.bnd;
+normals = bnd_mesh.nrm;
+centers = bnd_mesh.ctr;
+title('Boundary mesh');
+plot(bnd_mesh);
+hold on;
+
+%% Creating projection matrix linking boundary triangles to boundary edges
+
+% Number of elements in boundary and volume mesh
+Nelt_b = size(bnd_mesh.elt,1);
+Nelt_v = size(mesh.elt,1);
+tr_opr = sparse(Nelt_b, Nelt_v);
+
+% Looping over all boundary elements and finding the parent volume element
+for i = 1 : Nelt_b
+    bvindx = bnd_mesh.elt(i,:);
+    % Nodes of the boundary element
+    bvpts = bnd_mesh.vtx(bvindx,:);
+    % Checking which volume mesh nodes match with selected boundary nodes
+    distpt1 = vecnorm(mesh.vtx-ones(size(mesh.vtx,1),1).*bvpts(1,:),2,2);
+    distpt2 = vecnorm(mesh.vtx-ones(size(mesh.vtx,1),1).*bvpts(2,:),2,2);
+    % Indices for matching nodes
+    j1 = find(distpt1 == 0);
+    j2 = find(distpt2 == 0);
+    % Checking which volume element contains both nodes
+    for j = 1:Nelt_v
+       cur_elt = mesh.elt(j,:);
+       % Checking if coincides
+       check1 = cur_elt(1) == j1 | cur_elt(2) == j1 | cur_elt(3) == j1;
+       check2 = cur_elt(1) == j2 | cur_elt(2) ==j2 | cur_elt(3) == j2;
+       if check1 & check2
+           disp('Found match');
+          tr_opr(i,j)=1; 
+       end
+    end
+end
+quiver(centers(:,1),centers(:,2),normals(:,1),normals(:,2));
+
+%% Setting up the discrete system
+
+% Quadrature order
+order = 3;
+
+% Creating the integration domain
+Omega = dom(mesh, order);
+Gamma = dom(bnd_mesh, order);
+
+% FEM space without the dirichlet boundary condition
+S1_Omega = fem(mesh,'P1');
+% FEM space with the homogeneous Dirichlet BC
+S10_Omega = dirichlet(S1_Omega, bnd_mesh);
+S1_Gamma = fem(bnd_mesh,'P1');
+
+% LHS matrix
+A = integral(Omega,grad(S10_Omega),grad(S10_Omega));
+
+% Restriction operator for volume dofs to boundary dofs
+P = restriction(S1_Omega,bnd_mesh); % Gives a matrix with size S1_Gamma.Ndof X S1_Omega.Ndof
+
+% Getting the RHS vector
+g_N = g(S1_Gamma.dof);
+G_N = P' * g_N; % Extension by zero
+
+v = -integral(Omega,grad(S10_Omega),grad(S1_Omega))* G_N;
+
+% This is the offset solution u* 
+sol = A\v; 
+% Sol contains N_space - N_bnd elements, how to get a solution of size N_space?
+% Using Martin's approach of "elimination"
+Q = elimination(S1_Omega,bnd_mesh); % Matrix of size S1_Omega.Ndof X (S1_Omega.Ndof-S1_Gamma.Ndof)
+
+% Solution with full size N_space, extension by zero
+sol_ex = Q * sol;
+
+% Obtaining the solution u by undoing the offset
+u = sol_ex + G_N;
+
+%% TODO: Computing the solution to the adjoint equation
+
+
+%% Computing T1,T2,T3,...
+T1 = integral(Gamma,
+%%
+
+
+%% Plotting the state solution
+pts = S1_Omega.dof;
+X = pts(:,1);
+Y = pts(:,2);
+figure;
+scatter3(X,Y,u);
+figure;
+r = vecnorm(pts,2,2);
+u_exact = log(r/1.5)/log(0.5/1.5);
+scatter3(X,Y,abs(u_exact-u));
+title('Solution error');