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

addpath 'D:\EM_work\MT3D\mesh\Modeling script\lib'  ;
outname  = "ModelingTest1";% Model output file name
topofile = 'topoxyz.txt';
data_dir = pwd ;

% Use blocks to define out-of-model frames

%Size
blk.Lx = [50000,10000]; % unit：m
blk.Ly = [50000,10000];
blk.Lz = [50000,10000]; 
%Center position coordinate
blk.CentCoord = [0,0,0;0,0,0];
%Contrarotate
blk.angle = 0; %unit：deg

%
% Load the coordinates of the receiving point

load('rec.mat')
lengthcurve = 100; %The width of the surface should be less than the distance from the measuring point
heightcurve = 500; %Surface height, must be greater than the interpolation function error

% Load Topographic elevation data

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

%
%% 2 Modeling and
%
%% 2.1 Start COMSOL Multiphysics  with MATLAB
% If COMSOL Multiphysics with MATLAB is started, comment the following line
% 如果COMSOL Multiphysics  with MATLAB已经启动,需要注释下面这行代码

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：Building the Model Framework. 建立模型框架
% AddRecParametricSurface：The auxiliary parametric surface for constructing the measurement points is added. 加入构建测点的辅助参数曲面

[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();

%
%% 3 generating meshs
%% 3.1 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();

%
%% 3.2 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);