-
Notifications
You must be signed in to change notification settings - Fork 2
/
bem.m
58 lines (41 loc) · 1.4 KB
/
bem.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
% This implements a simple, indirect boundary element method for the case of circle and L-Shape.
% This code was written for OCTAVE.
% To run in MATLAB you mast replace the lines involving integration in MatAssembly.m and Poteval.m
%
% PLEASE, DO NOT REDISTRIBUTE THIS CODE.
% I CANNOT PROVIDE SUPPORT FOR THIS CODE.
% I CANNOT GUARANTEE THAT THIS CODE WORKS.
%
% Written by Felix Wolf
clear all;
close all;
% The approximate amount of elements to start with. The real number
% will be determined by mkGeom.
num = 4;
lmax = 5;
% Choose Method: 'c' for collocation, 'g' for Galerkin
method = 'c'
% Choose Geometry and corresponding rhs. 'l' for LShape, 's' for circle
geometry = 'l';
% Initialize arrays for error and number of elements
err = [];
nums = [];
% Loop for mesh refinement
for l = 1 : lmax
% Makes the geometry
Geom = mkGeom(num,geometry);
% Update the number for the next iteration
num = num*2;
% Store the number of elements
nums = [nums,max(size(Geom))-1];
% Assemble the matrix
A = mkMat(Geom,method);
% Assemble the RHS
U = mkRHS(Geom,geometry,method);
% Direct solver (CG wont work for collocation)
X = A\U';
% Evaluates the representation Formula and finds maximal potential error
[ptX,ptX1,ptX2,ptY,ptY1,ptY2,val,val1,val2,err] = repformeval(err,geometry,Geom,X);
err
end
visualize(Geom,nums,geometry,ptX,ptX1,ptX2,ptY,ptY1,ptY2,val,val1,val2,err);