-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathforward_solver.m
49 lines (34 loc) · 1.32 KB
/
forward_solver.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
%The main idea is to solve the forward proble two times. One for normalization and one for calculations
function fem=forward_solver(input,mesh)
fem=[];
input.jacobian_flag=1;
% Calculate k and g
[mesh.k,mesh.g]=rwaven(mesh.probe_spacing,mesh.probe_spacing*mesh.max_n,4);
% meah.mean_res=10;
mesh.prop(:)=mesh.mean_res;
tic
fem=mes_control_fast(0,input,mesh,fem,0);
toc
%update now res_param
mesh.res_param1=input.bgr_res_param;
for i=1:mesh.num_param
for j=1:mesh.num_elements
if(mesh.icon(4,j)==abs(i)) ;mesh.prop(j)=mesh.res_param1(i); end
end
end
%solve again for the user defined resistivity
input.jacobian_flag=1;
tic
fem=mes_control_fast(2,input,mesh,fem,0);
toc
if input.sip_flag==0
final=[input.ax input.az input.bx input.bz input.mx input.mz input.nx input.nz real(fem.array_model_data)];
model=[mesh.param_x mesh.param_y real(mesh.res_param1)];
else
final=[input.ax input.az input.bx input.bz input.mx input.mz input.nx input.nz real(fem.array_model_data) imag(fem.array_model_data)];
model=[mesh.param_x mesh.param_y real(mesh.res_param1) imag(mesh.res_param1)];
end
save('forward.d','final','-ascii');
% final=[input.ax;input.ay;input.bx;input.by;input.mx;input.my;input.nx;input.ny;fem.array_model_data];
% save('model4.mod','model','-ascii');
end