function [ r, Times, OutputRs, Wellpara, trun  ] = data_standard(  )
tic;
%% 模型选取
% 1--2014, Monifar    2--2018, steady, Rao     3--2018, transient, Rao
% 4-- PEDFM, Tene, Jiang, Younis
modelflag=1;
%前处理综合
%% 基质网格参数输入并进行网格参数计算
% dx=[10*ones(1,100) ];nx=size(dx,2);
% dy=[10*ones(1,50) ];ny=size(dy,2);
% dz=[10*ones(1,4)];nz=size(dz,2);
% dx=[10*ones(1,50) ];nx=size(dx,2);
% dy=[10*ones(1,30) ];ny=size(dy,2);
% dz=[10*ones(1,3)];nz=size(dz,2);
dx=[10*ones(1,100) ];nx=size(dx,2);
dy=[10*ones(1,50) ];ny=size(dy,2);
dz=[10*ones(1,1)];nz=size(dz,2);
%% 裂缝参数输入
% 两种输入方式，1表示向量输入，2表示基准点、倾角、方位角、抬升角输入
input_style=2;
if input_style==1
% ==================向量输入====================================== %
%首先对矩形缝进行研究，裂缝参数包括三个向量和两个参数的取值范围，如果是其它类型缝，则是两个参数间的函数关系
%因此一个5行三列的矩阵可以确定一条裂缝
%向量分量可为非整数，以此保证参数范围的取值为整数即可，
%故可先随意确定为整数的参数范围，再根据缝长缝高去确定向量分量的取值
 f=[
     0,1200,0;10,-20,0;0,0,10;0,60,0;0,3,0;
     0,1200,0;10,20,0;0,0,10;0,15,0;0,3,0;
     150,1500,0;20,-10,0;0,0,10;0,20,0;0,3,0;
     550,1300,0;10,10,0;0,0,10;0,10,0;0,3,0;
     650,1400,0;10,-20,0;0,0,10;0,55,0;0,3,0;
     605,705,0;10,0,0;0,0,10;-20,20,0;0,3,0;

%  385,485,2;5,0,0;0,0,1;-6,12,0;0,24,0;
%     385,365,2;5,0,0;0,0,1;-12,20,0;0,24,0;
%     385,245,2;5,0,0;0,0,1;-10,20,0;0,24,0;
%      385,925,2;5,0,0;0,0,1;-9,18,0;0,20,0;
%      385,805,2;5,0,0;0,0,1;-8,14,0;0,20,0;
%      385,685,2;5,0,0;0,0,1;-12,20,0;0,24,0;
%      385,565,2;5,0,0;0,0,1;-10,20,0;0,22,0;
%      385,445,2;5,0,0;0,0,1;-6,12,0;0,24,0;
%      955,1305,5;5,0,0;0,0,1;-5,4,0;0,15,0;
%      955,1225,5;5,0,0;0,0,1;-7,5,0;0,15,0;
%      955,1145,5;5,0,0;0,0,1;-9,6,0;0,15,0;
%      955,1065,5;5,0,0;0,0,1;-7,6,0;0,15,0;
%      955,985,5;5,0,0;0,0,1;-6,4,0;0,15,0;
%      955,755,5;5,0,0;0,0,1;-5,4,0;0,15,0;
%      955,675,5;5,0,0;0,0,1;-8,6,0;0,15,0;
%      955,595,5;5,0,0;0,0,1;-8,6,0;0,15,0;
%      955,515,5;5,0,0;0,0,1;-10,7,0;0,15,0;
%      955,435,5;5,0,0;0,0,1;-7,6,0;0,15,0;
%      955,355,5;5,0,0;0,0,1;-9,6,0;0,15,0;
%      1635,1245,5;5,0,0;0,0,1;-10,2,0;0,20,0;
%      1635,1165,5;5,0,0;0,0,1;-16,3,0;0,20,0;
%      1635,1085,5;5,0,0;0,0,1;-15,3,0;0,20,0;
%      1635,1005,5;5,0,0;0,0,1;-12,2,0;0,20,0;
%      1635,925,5;5,0,0;0,0,1;-10,2,0;0,20,0;
%      1635,695,5;5,0,0;0,0,1;-10,2,0;0,20,0;
%      1635,615,5;5,0,0;0,0,1;-12,2,0;0,20,0;
 ];
% 椭圆缝
% fellip=[320,210,10;10,2,0;0,0,2;0,0,8;0,0,5;
%            280,300,10;10,2,0;0,0,2;0,0,10;0,0,4;];
%     300,205,10;10,0,0;0,0,2;0,0,10;0,0,5];
fellip=[];
else
% ==================工程应用输入====================================== %
%仅能刻画具有双对称性质的矩形缝或椭圆缝
% 基准点坐标（1），方位角（2），倾角（3），抬升角（4），缝长（椭圆长轴长）（5），缝高（椭圆短轴长）（6），类型1是矩形缝、2是椭圆缝（7）
    input_content={
                     [255,255,5],90,90,0,200,10,1;
                     [305,255,5],90,90,0,200,10,1;
                     [355,255,5],90,90,0,200,10,1;
                     [405,255,5],90,90,0,200,10,1;
                     [455,255,5],90,90,0,200,10,1;
                     [505,255,5],90,90,0,200,10,1;
                     [555,255,5],90,90,0,200,10,1;
                     [605,255,5],90,90,0,200,10,1;
                     [655,255,5],90,90,0,200,10,1;
                     [705,255,5],90,90,0,200,10,1;
                     [755,255,5],90,90,0,200,10,1;
                   };
[f,fellip]=sort_fracture(input_content);
end

% ==================裂缝情况整理====================================== %
m=size(f,1);
n=size(fellip,1);
nf=(m+n)/5;%裂缝条数
%% 渗流介质相关参数
% =====================基质孔渗====================================== %
% kx = 0.1*1e-3 * ones(nx*ny*nz, 1);%生成基质网格绝对渗透率x方向矩阵
% ky =0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
% kz = 0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
pori = 0.3 * ones(nx*ny*nz, 1);%生成基质网格孔隙度矩阵
prpor = 16;%基质孔隙度基准压力
cpor = 1.07e-4;%孔隙体积压缩系数
NTG= 1 * ones(nx*ny*nz, 1);
% SRV区域定义
base_x=23:78; nx_SRV=length(base_x);
base_y=15:36; ny_SRV=length(base_y);
base_z=1:1;nz_SRV=length(base_z);
SRV=[];
for k=1:nz_SRV
    for j=1:ny_SRV
        for i=1:nx_SRV
            SRV=[SRV;(base_z(k)-1)*nx*ny+(base_y(j)-1)*nx+base_x(i)];
        end
    end
end  
% 非SRV渗透率赋值
kx = 0.1*1e-3 * ones(nx*ny*nz, 1);%生成基质网格绝对渗透率x方向矩阵
ky =0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
kz = 0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
% SRV渗透率赋值
kx(SRV)=2*1e-3 * ones(length(SRV), 1);
ky(SRV)=2*1e-3 * ones(length(SRV), 1);
kz(SRV)=2*1e-3 * ones(length(SRV), 1);
% load ('NTG.mat');
% NTG=sort_input(NTG,nx,ny,nz);
% =====================裂缝孔渗====================================== %
Kf = 20000*1e-3*ones(1,nf);%裂缝渗透率
% Kf = [0*1e-3*ones(1,5) 10000*1e-3*ones(1,1);];%裂缝渗透率
Wf =1e-2*ones(1,nf);%裂缝宽度
Porf = 0.30*ones(1,nf);%裂缝孔隙度
prporf = 20;%裂缝孔隙度基准压力
cporf = 1.07e-4;%裂缝孔隙体积压缩系数

Rpt = 2;%flag
cf = 1;
ca = 1;
% co = 3.02e-3;
%% 流体性质参数
% =======================密度================================= %
Dosi=820; %Dosi是地面标况下测得的油相密度，单位是千克/每立方米
Dwsi=1000; %Dwsi是地面标况下测得的水相密度
% ====================体积系数及压缩系数================================= %
pb = 20.0;%泡点压力
Bopb = 1.0;%泡点压力pb对应的原油体积系数
co = 3.02e-3;%原油体积压缩系数
prw = 20.0;%地层水体积系数基准压力
Bwi = 1.00001;%原始地层水体积系数
cw = 5e-4;%地层水体积压缩系数
prg = 25.0;%地层水体积系数基准压力
% =======================粘度================================= %
visopb = 2;%泡点压力pb对应的原有粘度
cvo = 0;%原油粘度随压力变化的系数
vwi =0.6;%原始地层水粘度
cvw = 0;%地层水粘度随压力变化的系数
% ==================相渗及毛管力曲线====================================== %
% ==================基质油水相渗====================================== %
RPOW=[
0	0	0.8	0
0.15	0	0.8	0
0.45	0.04	0.4	0
0.5212	0.08	0.227	0
0.5425	0.103	0.184	0
0.5638	0.125	0.156	0
0.5851	0.145	0.125	0
0.6064	0.174	0.094	0
0.6276	0.207	0.063	0
0.6489	0.247	0.046	0
0.6701	0.292	0.031	0
0.6914	0.335	0.021	0
0.7127	0.378	0.011	0
0.76	0.475	0.005	0
0.85	0.552	0	0
1	0.552	0	0
];
SW=RPOW(:,1);KRW=RPOW(:,2);KRO=RPOW(:,3);PCOW=RPOW(:,4);
% ==================裂缝油水相渗====================================== %
% 在裂缝中，性质就是简单得两点线性插值
if Rpt ~= 1%表示该网格是裂缝网格
PRF=[
0	0	0.8	0
0.15	0	0.8	0
0.85	0.6	0	0
1	0.6	0	0
];
SWF = PRF(:,1);KRWF = PRF(:,2);KROF = PRF(:,3);PCOWF=PRF(:,4);
end
% ==================毛管力情况====================================== %
ifpcow = 0;%此处为控制是否忽略油水相间毛管力，0代表忽略，非0代表考虑
%% 结合网格和流动介质参数进行前处理
r = GridProp(modelflag, nx, ny, nz, dx, dy, dz, kx, ky, kz, f, fellip, Kf, Wf, pori, prpor, cpor,  Porf, prporf, cporf, cf, ca,co, NTG);
%% 生成向量化变成所需要的算子
os = OperatorRS(r.N, r.nex, r.nc);
%% 结合流体性质参数进行处理
f = fluidPVT(Bopb, pb, co, Bwi, prw, cw, vwi, cvw, visopb, cvo, SW, KRO, KRW, PCOW, ifpcow, SWF, KROF, KRWF, PCOWF, r.rpt,Dosi,Dwsi);
%% 初值条件 Initial Condition Section
% =================压力、饱和度初值======================================= %
%压力初值要考虑重力，给出油藏下表面压力值
P = [20 * ones(r.nmc, 1); 20 * ones(r.nfc, 1)];
% plow=25; % P =plow-1e-6*800*9.8*r.z;
Sw = [0.15 * ones(r.nmc, 1); 0.15 * ones(r.nfc, 1)];
state0 = initialRS(P, Sw);
%%  schedule
% 井参数输入
% Wellcpara: welltype(1) nperf(2) index(3) rw(4) protype(5) value(6)
% constraints(7) skin(8) welltype2(9) 
% welltype = 1 prod, welltype = 2 inj
% protype = 1 const flowrate, protype = 2 const pwf
% welltype2=1 直井， welltype2=2 水平井沿x方向，welltype2=3 水平井沿y方向 welltype2=4 多段压裂水平井 
%将直井、水平井与多段压裂水平井分开处理
%直井、水平井，该种井型的处理是将射孔段安排在基质网格
%由于是三维的储层，具体基质网格不好直接给出，对于直井、斜井、水平井，采取与eclipse一致的方式，给出其在基质网格中位置
%每个射孔点由下px,py,pz确定，分别表示该射孔点在x方向，y方向的网格编号，及层数即z反方向上的网格编号（最初是以垂直向上建立的基质网格）
well1={
%            2, 3, [10 20 1;10 20 2;10 20 3], 0.178/2, 2, 35 , 35, 0, 1;
%            2, 3, [67 20 1;67 20 2;67 20 3], 0.178/2, 2, 35 , 35, 0, 1;
%            2, 3, [130 20 1;130 20 2;130 20 3], 0.178/2, 2, 35 , 35, 0, 1;
%            2, 3, [196 20 1;196 20 2;196 20 3], 0.178/2, 2, 35 , 35, 0, 1;
%            2, 3, [10 148 1;10 148 2;10 148 3], 0.178/2, 2, 35 , 35, 0, 1;
%            2, 3, [67 148 1;67 148 2;67 148 3], 0.178/2, 2, 35 , 35, 0, 1;
%            2, 3, [130 148 1;130 148 2;130 148 3], 0.178/2, 2, 35 , 35, 0, 1;
%            2, 3, [196 148 1;196 148 2;196 148 3], 0.178/2, 2, 35 , 35, 0, 1;
           };
%需要找到射孔段具体所在的基质网格编号
n1=size(well1,1);
for i=1:n1
    perf=well1{i,3};
    nperf=well1{i,2};
    perf(:,3)=(r.nz+1)*ones(nperf,1)-perf(:,3);
    A=repmat([0 1 1],nperf,1);
    B=(perf-A)*[1;r.nx;r.nx*r.ny];
    B=B';    %转换为行向量
    well1{i,3}=B;
end
% %多段压裂水平井的处理，射孔段设置在裂缝单元上，可以给出该射孔点所在的坐标，然后去寻找该点所在的裂缝单元编号
welloc1=[
    [255,255,5];
                     [305,255,5];
                     [355,255,5];
                     [405,255,5];
                     [455,255,5];
                     [505,255,5];
                     [555,255,5];
                     [605,255,5];
                     [655,255,5];
                     [705,255,5];
                     [755,255,5];
%     [205,255,15];
%                    [405,255,21];
%                    [605,255,21];
%                    [805,255,31];
%              385,805,15;
%              385,685,15;
%              385,565,15;
%              385,445,15;
             ];
% welloc2=[955,1305,15;
%                955,1225,15;
%                955,1145,15;
%                955,1065,15;
%                955,985,15;
%                955,755,15;
%                955,675,15;
%                955,595,15;
%                955,515,15;
%                955,435,15;
%                955,355,15;];
% welloc3=[1635,1245,15;
%          1635,1165,15;
%          1635,1085,15;
%          1635,1005,15;
%          1635,925,15;
%          1635,695,15;
%          1635,615,15;
%                ];
perfnum1 = findWelloc(r, welloc1);%在用以Wellc0之中
% perfnum2 = findWelloc(r, welloc2);
% perfnum3 = findWelloc(r, welloc3);
well2={
    1, length(perfnum1), perfnum1, 0.178/2, 2,10, 10, 0,4;
%     1, length(perfnum2), perfnum2, 0.178/2, 2,8, 8, 0,4;
%     1, length(perfnum3), perfnum3, 0.178/2, 2,8, 8, 0,4;
           };
% well2={};
% well2={1, length(perfnum), perfnum, 0.1, 2,10, 10, 0,4;
%            };
%        2, 1, [  5030 ], 0.3/0.318, 1,0.05, 30, 0,4
%把两类井整合到一起
nwell1=size(well1,1);
nwell2=size(well2,1);
nwell=nwell1+nwell2;
Wellc0=cell(nwell,9);
if nwell1==0
    Wellc0=well2;
elseif nwell2==0
    Wellc0=well1;
else
for i=1:nwell1
    for j=1:9
    Wellc0{i,j}=well1{i,j};
    end
end
for i=1:nwell2
    for j=1:9
    Wellc0{i+nwell1,j}=well2{i,j};
    end
end
end
 
yitap = 5; % 50-500psi
yitas = 0.04; % 0.05-0.5
omega = 0.5; % 0-1
Nmax = 50;
epsave = 1e-6;
epsmax = 1e-6;
dtmin = 0.001;
dtmax =20;
tend = 500;
nwel = size(Wellc0, 1);
WelChg = zeros(nwel, 1);%这表示都开井
Wellc = calcTrans(r, Wellc0);
w.yitap = yitap;
w.yitas = yitas;
w.omega = omega;
w.Nmax = Nmax;
w.epsave = epsave;
w.epsmax = epsmax;
w.epsave = epsave;
w.epsmax = epsmax;
w.dtmin = dtmin;
w.dtmax = dtmax;
w.tend = tend;
w.WelChg = WelChg;
w.Wellc = Wellc; 

%% 计算
if modelflag==1 || modelflag==4
    [Times, OutputRs, Wellpara, trun] = mainRS_MB_2014(r, f, os, w, state0);
else
[Times, OutputRs, Wellpara, trun] = mainRS_MB_modified(r, f, os, w, state0);
end
% [Times, OutputRs, Wellpara, trun] = mainRS_MB_modified(r, f, os, w, state0);
run_time=toc;
trun.run_time=run_time;
  end