code

Signed-off-by: 刘明宏 <liuminghong@mail.sdu.edu.cn>

bin/RMT3d/GetInterpMatrix.m

@@ -0,0 +1,85 @@
+function Q = GetInterpMatrix(rec,mesh,indelem2)
+%% Field interpolation function
+
+% rec  = [0 0 0.20];
+% rec  = randn(1000,3);
+
+t = tic;
+mu0 = 4*pi*1e-7  ;  % V·s/ (A·m)
+nrec = size(rec,1);
+if nargin==2
+    TR = triangulation(mesh.elem2node,mesh.node2coord);
+    ID = pointLocation(TR,rec);
+elseif nargin==3
+    ID = indelem2;
+end
+
+a           = mesh.abcd(ID,[1   2  3  4]);
+b           = mesh.abcd(ID,[5   6  7  8]);
+c           = mesh.abcd(ID,[9  10 11 12]);
+d           = mesh.abcd(ID,[13 14 15 16]);
+ve          = mesh.vol(ID,1);
+Sign        = mesh.sign(ID,:);
+elem2edge   = mesh.elem2edge(ID,:);
+elemEdgeLength  = mesh.elemEdgeLength(ID,:);
+local_i    = [1 2;1 3;1 4;2 3;4 2;3 4];
+
+[Nie, CurlNie] = getEit1(rec(:,1),rec(:,2),rec(:,3),a,b,c,d,ve,local_i,Sign,elemEdgeLength);
+
+
+Q.Ex=sparse(nrec,mesh.nedge);
+Q.Ey=sparse(nrec,mesh.nedge);
+Q.Ez=sparse(nrec,mesh.nedge);
+Q.Hx=sparse(nrec,mesh.nedge);
+Q.Hy=sparse(nrec,mesh.nedge);
+Q.Hz=sparse(nrec,mesh.nedge);
+
+ii=repmat((1:nrec)',1,6);
+
+Q.Ex=sparse(ii,elem2edge,squeeze(Nie(:,1,:)),nrec,mesh.nedge);
+Q.Ey=sparse(ii,elem2edge,squeeze(Nie(:,2,:)),nrec,mesh.nedge);
+Q.Ez=sparse(ii,elem2edge,squeeze(Nie(:,3,:)),nrec,mesh.nedge);
+Q.Hx=sparse(ii,elem2edge,squeeze(CurlNie(:,1,:)),nrec,mesh.nedge);
+Q.Hy=sparse(ii,elem2edge,squeeze(CurlNie(:,2,:)),nrec,mesh.nedge);
+Q.Hz=sparse(ii,elem2edge,squeeze(CurlNie(:,3,:)),nrec,mesh.nedge);
+
+Q.Hx = -Q.Hx;
+Q.Hy = -Q.Hy;
+Q.Hz = -Q.Hz;
+
+Q.Hx=Q.Hx/mu0;
+Q.Hy=Q.Hy/mu0;
+Q.Hz=Q.Hz/mu0;
+
+GetInterpMatrixTime = toc(t)
+end
+
+
+
+function [Nie, CurlNie]=getEit1(xc,yc,zc,a,b,c,d,ve,local_i,TE_Sign,TE_length)
+Ne          = size(xc,1);
+Nie         = zeros(Ne,3,6);
+CurlNie     = zeros(Ne,3,6);
+
+for kt1 = 1:6
+    i1              = local_i(kt1,1);
+    i2              = local_i(kt1,2);
+    Sign_i          = TE_Sign(:,kt1);
+    ai1             = a(:,i1);
+    ai2             = a(:,i2);
+    bi1             = b(:,i1);
+    bi2             = b(:,i2);
+    ci1             = c(:,i1);
+    ci2             = c(:,i2);
+    di1             = d(:,i1);
+    di2             = d(:,i2);
+    l               = TE_length(:,kt1);
+    Li1             = (ai1 + bi1.*xc + ci1.*yc + di1.*zc)./(6*ve);
+    Li2             = (ai2 + bi2.*xc + ci2.*yc + di2.*zc)./(6*ve);
+    GradLi1         = [bi1,ci1,di1]./repmat(6*ve,1,3);
+    GradLi2         = [bi2,ci2,di2]./repmat(6*ve,1,3);
+    Nie(:,:,kt1)    = repmat(Sign_i,1,3).*(repmat(Li1,1,3).*GradLi2 - repmat(Li2,1,3).*GradLi1).*repmat(l,1,3);
+    CurlNie(:,:,kt1) =  repmat(Sign_i,1,3).*[ci1.*di2-di1.*ci2, di1.*bi2-bi1.*di2, bi1.*ci2-ci1.*bi2].*[2*l./((6*ve).^2),2*l./((6*ve).^2),2*l./((6*ve).^2)];
+end
+
+end

bin/RMT3d/Set_material_3Ani.m

@@ -0,0 +1,35 @@
+function [parameter] = Set_material_3Ani(entity,anoID,airID,ano_xx,ano_yy,ano_zz,air,parameter_residual)
+% entity: nelm*1, the domain index of a tetrahedral unit
+% anoID: the domain index of anomalous body
+% airID: the domain index of air
+% ano_xx,ano_yy,ano_zz: the anisotropy parameter of anomalous body
+% air:  parameter of air domain
+% parameter_residual: parameter of residual domain
+nelems     = size(entity,1);
+parameter_xx   = zeros(nelems,1);       
+parameter_yy   = zeros(nelems,1);
+parameter_zz   = zeros(nelems,1);
+na      = length(anoID);    %Number of anomalous domain
+n_airentity=length(airID);
+
+for i=1:n_airentity
+    %设置空气物性参数 Set the physical property parameters of air
+    parameter_xx(entity == airID(i))= air;
+    parameter_yy(entity == airID(i))= air;
+    parameter_zz(entity == airID(i))= air;
+end
+
+% 设置anoID区域物性参数 Set the physical property parameters of anoID
+for ia = 1:na
+    parameter_xx(entity==anoID(ia))=ano_xx(ia);
+    parameter_yy(entity==anoID(ia))=ano_yy(ia);
+    parameter_zz(entity==anoID(ia))=ano_zz(ia);
+end
+
+ %设置剩余区域物性参数 Set the physical property parameters of residual domain
+parameter_xx(~ismember(entity, [anoID airID]))= parameter_residual;
+parameter_yy(~ismember(entity, [anoID airID]))= parameter_residual;
+parameter_zz(~ismember(entity, [anoID airID]))= parameter_residual;
+
+parameter = [parameter_xx,parameter_yy,parameter_zz];
+end

bin/RMT3d/get_MTR.m

@@ -0,0 +1,55 @@
+function [MTR] = get_MTR(Ex1,Ey1,Hx1,Hy1,Hz1,Ex2,Ey2,Hx2,Hy2,Hz2,freq)
+
+% Magnetotelluric response calculation
+% Return: impedance, apparent resistivity, phase,  tipper response
+
+mu0 = 4*pi*1e-7 ;  % V·s/ (A·m)
+
+
+omega = 2*pi*freq;
+detH  = (Hx1.*Hy2-Hx2.*Hy1);
+
+%% impedance
+Zxx  = (Ex1.*Hy2 - Ex2.*Hy1)./detH ;
+Zxy  = (Ex2.*Hx1 - Ex1.*Hx2)./detH ;
+Zyx  = (Ey1.*Hy2 - Ey2.*Hy1)./detH ;
+Zyy  = (Ey2.*Hx1 - Ey1.*Hx2)./detH ;
+
+MTR.Zxx = Zxx;
+MTR.Zxy = Zxy;
+MTR.Zyx = Zyx;
+MTR.Zyy = Zyy;
+
+%% apparent resistivity
+rhoXX  =  1/(omega*mu0)*abs(Zxx).^2;
+rhoXY  =  1/(omega*mu0)*abs(Zxy).^2;
+rhoYX  =  1/(omega*mu0)*abs(Zyx).^2;
+rhoYY  =  1/(omega*mu0)*abs(Zyy).^2;
+
+MTR.rhoXX = rhoXX;
+MTR.rhoXY = rhoXY;
+MTR.rhoYX = rhoYX;
+MTR.rhoYY = rhoYY;
+%% phase
+phaXX   =  -atan(imag(Zxx)./real(Zxx))*180/pi;
+phaXY   =  -atan(imag(Zxy)./real(Zxy))*180/pi;
+phaYX   =  -atan(imag(Zyx)./real(Zyx))*180/pi;
+phaYY   =  -atan(imag(Zyy)./real(Zyy))*180/pi;
+
+% phaXX   =  -atan2(imag(Zxx),real(Zxx))*180/pi;
+% phaXY   =  -atan2(imag(Zxy),real(Zxy))*180/pi;
+% phaYX   =  -atan2(imag(Zyx),real(Zyx))*180/pi;
+% phaYY   =  -atan2(imag(Zyy),real(Zyy))*180/pi;
+MTR.phaXX = phaXX;
+MTR.phaXY = phaXY;
+MTR.phaYX = phaYX;
+MTR.phaYY = phaYY;
+
+%% Tipper
+Tzx       = (Hz1.*Hy2 - Hz2.*Hy1)./detH ;
+Tzy       = (Hz2.*Hx1 - Hz1.*Hx2)./detH ;
+MTR.Tzx = Tzx;
+MTR.Tzy = Tzy;
+
+end
+

bin/RMT3d/load_comsol.m

@@ -0,0 +1,6 @@
+function [p,t,entity,sort_ind_sign] = load_comsol(name)
+load(name);
+p = comsolpt.p;
+t = comsolpt.t;
+entity = comsolpt.entity;
+% sort_ind_sign= comsolpt.sort_ind_sign;%场源棱边编号 CSEM、TEM

bin/RMT3d/load_mesh.m

@@ -0,0 +1,4 @@
+
+mesh.entity =load('mesh.entity.dat');
+mesh.elem2node =load('mesh.elem2node.dat');
+mesh.node2coord =load('mesh.node2coord.dat');

bin/RMT3d/save_mesh_mat2dat.m

@@ -0,0 +1,17 @@
+function  save_mesh_mat2dat(mesh)
+%save_mesh_mat2dat .mat to .dat
+%   此处显示详细说明
+
+entity = double(mesh.entity);  
+try elem2node = double(mesh.t); end 
+try elem2node = double(mesh.elem2node); end
+try node2coord = double(mesh.p); end
+try node2coord = double(mesh.node2coord); end
+
+save('mesh.entity.dat', 'entity', '-ascii');
+save('mesh.elem2node.dat', 'elem2node', '-ascii');
+save('mesh.node2coord.dat', 'node2coord', '-ascii');
+
+% mesh.entity =load('mesh.entity.dat');
+% mesh.elem2node =load('mesh.elem2node.dat');
+% mesh.node2coord =load('mesh.node2coord.dat');

bin/RMT3d/signs_edges.m

@@ -0,0 +1,43 @@
+function signs = signs_edges(elems2nodes)
+
+% SIGNS_EDGES
+%   Calculate signs for the edges of a 2D/3D mesh.
+%
+%   This data is needed in order to use linear Nedelec's elements in 2D/3D
+%   and the linear Raviart-Thomas element in 2D. For the linear Raviart-
+%   Thomas element in 3D we need signs related to faces, and this data is
+%   provided by the function 'signs_faces(...)'.
+%
+%   The edge signs can be easily deduced from the mesh data itself by
+%   directly using the data structure which represents the elements
+%   (triangles or tetrahedrons) by their nodes. The signs are obtained
+%   with minimal matrix operations in a vectorized manner.
+%
+% SYNTAX:  signs = signs_edges(elements)
+%
+% IN:   elems2nodes  Elements by their nodes.
+%                    In 2D, elements(i,1:3) are the three nodes of the i'th
+%                    triangle, and in 3D elements(i,1:6) are the six nodes of
+%                    the i'th tetrahedron.
+%
+% OUT:  signs        Signs for element edges, corresponding to the data
+%                    structure 'elements': signs(i,j) is the sign related to
+%                    the j'th edge of the i'th triangle (or tetrahedron).
+%
+
+dim = size(elems2nodes,2)-1;
+
+if ( dim == 2 )
+    tmp = elems2nodes(:,[2 3 1]) - elems2nodes(:,[3 1 2]);
+    signs = tmp ./ abs(tmp);
+elseif (dim == 3)
+    %     tmp = elems2nodes(:,[1 1 1 2 3 4]) - elems2nodes(:,[2 3 4 3 4 2]);
+    % **** 20221026
+    tmp = elems2nodes(:,[1 1 1 2 4 3]) - elems2nodes(:,[2 3 4 3 2 4]);
+
+    signs = tmp ./ abs(tmp);
+else
+    error('The data is not understood.')
+end
+
+end

bin/RMT3d/umfpack_MT.m

@@ -0,0 +1,26 @@
+function [x,res] = umfpack_MT(A,b)
+% Solve for x in linear system A * x = b.
+
+fprintf ('Solution to Ax=b via UMFPACK\n') ;
+
+% Compute the numeric factorization via UMFPACK.
+t = tic;
+[L,U,P,Q,R,Info] = umfpack (A, struct ('details',0));
+time(1) = toc(t);
+fprintf (['The LU factorization time is ' num2str(time(1)) 's\n']) ;
+
+% Compute the solution x using the symbolic factorization.
+t = tic;
+x = Q*(U\(L\(P*(R\b))));
+x = full(x);
+time(2) = toc(t);
+fprintf (['The symbolic factorization time is ' num2str(time(2)) 's\n']) ;
+
+% Compute the residuals.
+
+for i = 1:2%size(b,2)
+    res(i) = norm(A*x(:,i)-b(:,i))/norm(b(:,i));
+    fprintf('%s %d\n',['Normalized residual from SuiteSparse for rhs' num2str(i) ' is'], res(i));
+end
+
+end