%% GeoMesh
% Generate a tetrahedral mesh with topography, this interface is based on  COMSOL Multiphysics with MATLAB
%% 1 Model parameter definition
%
%% Model output file name

clear
addpath 'D:\EM_work\MT3D\mesh\Modeling script\lib'
outname  = 'KarsTunnel0';
topofile = 'topoxyz.txt';
data_dir  = pwd ;
RX       = 50 ;               % Receiving distance
p        = textread(topofile);
load rec.mat;

% Use blocks to define out-of-model frames

%Size
Lx1 = max(p(:,1))-min(p(:,1));
Lx2 = max(rec(:,1))-min(rec(:,1))+5*RX;
Ly1 = max(p(:,2))-min(p(:,2));
Ly2 = max(rec(:,2))-min(rec(:,2))+5*RX;
blk.Lx = [Lx1,Lx2]; % unit：m
blk.Ly = [Ly1,Ly2] ;
blk.Lz = [5000,1000] ; 
%Center position coordinate
recCen = mean(rec) ;
outCen = mean(p);
blk.CentCoord = [outCen;recCen];
%Contrarotate
blk.angle = 0; %unit：deg

%
%% Defined receiving point

lengthcurve = 0.25*RX; %The width of the surface should be less than the distance from the measuring point
heightcurve = RX;      %Surface height, must be greater than the interpolation function error

%% Defined topographic surface

minX = outCen(1)-0.5*blk.Lx(:,1);
maxX = outCen(1)+0.5*blk.Lx(:,1);
minY = outCen(2)-0.5*blk.Ly(:,1);
maxY = outCen(2)+0.5*blk.Ly(:,1);

%% Generative model geometry
%% 2.1 启动COMSOL Multiphysics 6.0 with MATLAB
% 如果COMSOL Multiphysics  with MATLAB已经启动,需要注释下面这行代码
% If COMSOL Multiphysics with MATLAB is started, comment the following line

 Comsol_with_Matlab_Start;

%% 2.2  Modeling
% 初始化 3D model initialization

import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');% ModelUtil.remove('Model')
model.modelNode.create('mod1');
model.geom.create('geom1', 3);
model.mesh.create('mesh1', 'geom1');

% AddBlocks建立模型框架
% AddRecParametricSurface加入测点参数曲面

[object_blk,blk_pos_corner ] = AddBlocks(model,blk) ;
[object_rec,ps_pos_rec ] = AddRecParametricSurface(model,rec,lengthcurve,heightcurve);

%
% Add topography parameter surface

model.component('mod1').func.create('int1', 'Interpolation');
model.component('mod1').func('int1').set('source', 'file');
model.component('mod1').func('int1').set('importedname', topofile);
model.component('mod1').func('int1').set('importedstruct', 'Spreadsheet');
model.component('mod1').func('int1').set('importeddim', '2D');
model.component('mod1').func('int1').set('argunit', {'' ''});
model.component('mod1').func('int1').set('filename', [data_dir '\' topofile]);
model.component('mod1').func('int1').importData;
model.component('mod1').func('int1').set('nargs', '2');
model.component('mod1').func('int1').set('struct', 'spreadsheet');

model.component('mod1').geom('geom1').create('pstopo', 'ParametricSurface');
model.component('mod1').geom('geom1').feature('pstopo').set('parmin1', minX);
model.component('mod1').geom('geom1').feature('pstopo').set('parmax1', maxX);
model.component('mod1').geom('geom1').feature('pstopo').set('parmin2', minY);
model.component('mod1').geom('geom1').feature('pstopo').set('parmax2', maxY);
model.component('mod1').geom('geom1').feature('pstopo').set('coord', {'s1' 's2' 'int1(s1,s2)'});
model.component('mod1').geom('geom1').feature('pstopo').set('maxknots', 100);
model.component("mod1").geom("geom1").run();

% Convert to Solid

object = [object_blk, object_rec,{'pstopo'}];
csol1  = model.component("mod1").geom("geom1").create("csol1", "ConvertToSolid");
csol1.selection('input').set(object);
model.component("mod1").geom("geom1").run();

%
% Remove the superfluous boundary introduced at the measuring point

idBound_all=[];
for nrec=1:size(ps_pos_rec,1)
idx1 = mphselectcoords(model,'geom1',ps_pos_rec(nrec,1:3),'boundary','radius',0.1);
idx2 = mphselectcoords(model,'geom1',ps_pos_rec(nrec,4:6),'boundary','radius',0.1);
idBound_all = [idBound_all,[idx1 idx2]];
end
del = model.component("mod1").geom("geom1").create("del1", "Delete");
del.selection('input').set('csol1',idBound_all);
model.component("mod1").geom("geom1").run();

% Lofting the shape of the cave

CurveFiles =  get_curveTxtFile(data_dir,[1:5],'curve');
objectIPC = AndInterpolationCurve(model,CurveFiles);

model.geom("geom1").create("loft1", "Loft");
model.geom("geom1").feature("loft1").selection("profile").set(objectIPC);
model.geom("geom1").feature("loft1").set("facepartitioning", "grid");
model.component("mod1").geom("geom1").run();

%
%% 2.3 generate tetrahedral mesh
% The default parameters are used to generate tetrahedral mesh

model.component("mod1").mesh("mesh1").create("ftet1", "FreeTet");
model.component("mod1").mesh("mesh1").run();

%
%% 2.4 Output mesh information

[stats,data]        =mphmeshstats(model, 'mesh1');
comsolpt.entity             =data.elementity{2};
entity             =data.elementity{2};
t                  = data.elem{2}'+1;
p                  = data.vertex';
p(:,end)           = -p(:,end);
t                  = double(t);

 comsolpt.p = p;
 comsolpt.t = t;
 comsolpt.entity = entity;
 save(outname,'comsolpt');

 mphsave(model,outname);