~wbetz/fesslix/flx_testing

# Copying and distribution of this file, with or without modification,
# are permitted in any medium without royalty provided the copyright
# notice and this notice are preserved.  This file is offered as-is,
# without any warranty.

# data_import - mesh&RF (SOFiSTiK)
# ================================

loadlib "rnd";
default log::level = 4;

const PO = 1;
const PPDIM_EOLE = 5;

# Define the random field (Eur)  -> log-normal
# -----------------------------

const SD_1 = 25600;
const MU_1 = 8e4;
const a_1 = 5;
var cor_1 = exp(-1*(deltap/a_1)^2);
  const v_1 = SD_1/MU_1;
  const R_1 = v_1/sqrt(ln(1+V_1^2));
  rnd rf new 1 kl_fem_Gauss mu 0 sd 1 rhocor cor_1*R_1;

data_import("test_25_mesh.dat",SOFiSTiK) 
  (rf_elset: 1 pOrderLx=PO pOrderLy=PO eole_smp_mth=1 eole_n_lx=PPDIM_EOLE/2 eole_n_ly=PPDIM_EOLE 
  INNER_SAME_MIN_LEVEL=1 INNER_DIFF_MIN_LEVEL=0 OUTER_MIN_LEVEL=0 ERR_ID_MIN_LEVEL=0 ERR_IS_MIN_LEVEL=1
  );

# Define the random field (Eod)  -> log-normal
# -----------------------------

const SD_2 = 9600;
const MU_2 = 3e4;
const a_2 = 5;
var cor_2 = exp(-1*(deltap/a_2)^2);
  const v_2 = SD_2/MU_2;
  const R_2 = v_2/sqrt(ln(1+V_2^2));
  rnd rf new 2 kl_fem_Gauss mu 0 sd 1 rhocor cor_2*R_2;

data_import("test_25_mesh.dat",SOFiSTiK) 
  (rf_elset: 2 pOrderLx=PO pOrderLy=PO eole_smp_mth=1 eole_n_lx=PPDIM_EOLE/2 eole_n_ly=PPDIM_EOLE 
  INNER_SAME_MIN_LEVEL=1 INNER_DIFF_MIN_LEVEL=0 OUTER_MIN_LEVEL=0 ERR_ID_MIN_LEVEL=0 ERR_IS_MIN_LEVEL=1
  );


const lambda_1 = ln(MU_1)-0.5*ln(SD_1^2/MU_1^2+1);
  const chi_1 = sqrt(ln(SD_1^2/MU_1^2+1));
const lambda_2 = ln(MU_2)-0.5*ln(SD_2^2/MU_2^2+1);
  const chi_2 = sqrt(ln(SD_2^2/MU_2^2+1));
calc exp(0*chi_1+lambda_1);
calc exp((0.7*0+sqrt(1-0.7^2)*0)*chi_2+lambda_2);


# Solve the RF-discretization
# ---------------------------

const DIM = leak_check?5:50;
  rnd rf solve 1 dim=DIM;
  rnd rf solve 2 dim=DIM;
  sfor (i;DIM) { 
    funplot i,
          {tmpmtx!{
            mtxconst_new tmpmtx = {1,2};
            mtxconst_op tmpmtx(j) = rfeigenvalue(j,i);
          }};
  };
if (leak_check) {} else {
 const PPDIM = 1;
  filestream fs_rf_p_1_e1 ( "output/test_25_rf_1_plot_e1.txt") ;
  rnd rf plot 1 { stream = fs_rf_p_1_e1; ppdim={PPDIM;PPDIM}; };
};


# Test the element-integral
# -------------------------

const NGP = 10;
qid = 20018;

var rlz_1 = exp(rf(1)*chi_1+lambda_1);
var rlz_2 = exp((0.7*rf(1)+sqrt(1-0.7^2)*rf(2))*chi_2+lambda_2);

var av_1 = geointeg(face,qid,rlz_1,true,NGP);
var av_2 = geointeg(face,qid,rlz_2,true,NGP);
calc av_1;
calc av_2;

rnd smp { dolog=false; verbose=false; };
calc av_1;
calc av_2;

rnd smp { dolog=false; verbose=false; };
calc av_1;
calc av_2;

rnd smp { dolog=false; verbose=false; };
calc av_1;
calc av_2;

rnd smp { dolog=false; verbose=false; };
calc av_1;
calc av_2;

rnd smp { dolog=false; verbose=false; };
calc av_1;
calc av_2;


# Prepare a SOFiSTiK-File
# -----------------------

filestream sstream ( "output/test_25_system_p2.dat" );
filestream mstream ( "output/test_25_mat_p1.dat" );

file_filter_sofistik("test_25_system.dat") sstream, mstream, qid, mid, " MAT  @{%mid} E @{av_1} MUE $(L2_MUE) GAM $(L2_GAM) TITL 'Soil: Layer 2'
  NMAT @{%mid} GRAN P1 $(L2_FRC) P2 $(L2_COH) P5 @{av_2} P9 @{av_2} P10 $(L2_m)", {2}, 10 {};

filestream sstream ( "output/test_25_system_p1.dat" );

file_filter_sofistik("output/test_25_system_p2.dat") sstream, mstream, qid, mid, " MAT  @{%mid} E @{av_1} MUE $(L2_MUE) GAM $(L2_GAM) TITL 'Soil: Layer 2'
  NMAT @{%mid} GRAN P1 $(L2_FRC) P2 $(L2_COH)+$(COHPNL) P5 @{av_2} P9 @{av_2} P10 $(L2_m)", {3}, mid+1 {};

filestream sstream ( "output/test_25_system.dat" );
ostream_close mstream;

file_filter_cv ("output/test_25_system_p1.dat") { stream = sstream; };

run_external "rm output/test_25_system.dat";
run_external "rm output/test_25_mat_p1.dat";


# 'Post-processing'
# -----------------------

const N_smp = 10;
const N_el = 923;

sfor (i; N_smp) {
  rnd smp { dolog=false; verbose=false; };
  filestream res ("output/test_25_tmp_" & {i}$int ) { dolog=false; };
  filter(qid;tmpvec!{ mtxconst_rf_elidvec tmpvec(1,face);}) {
    funplot av_1 { stream=res; };
  };
  filestream res2 ("output/test_25_tmp2_" & {i}$int ) { dolog=false; };
  filter(qid;tmpvec) {
    funplot av_2 { stream=res2; };
  };
};

# Get the mean
  mtxconst_new postexp(N_el,1,0);
  mtxconst_new tmpmtx2({postexp});
  sfor (i; N_smp) {
    ifstream irnd ( "output/test_25_tmp_" & {i}$int ) { dolog=false; };
    mtxconst_fromfile tmpmtx2(col=1)=irnd;
    mtxconst_op postexp += {tmpmtx2};
  };
  mtxconst_op postexp /= N_smp;
  calc mtxsum(postexp)/N_el;
  mtxconst_new postexp2(N_el,1,0);
  sfor (i; N_smp) {
    ifstream irnd ( "output/test_25_tmp2_" & {i}$int ) { dolog=false; };
    mtxconst_fromfile tmpmtx2(col=1)=irnd;
    mtxconst_op postexp2 += {tmpmtx2};
  };
  mtxconst_op postexp2 /= N_smp;
  calc mtxsum(postexp2)/N_el;
# Get the standard deviation
  mtxconst_new postsd(N_el,1,0);
  sfor (i; N_smp) {
    ifstream irnd ( "output/test_25_tmp_" & {i}$int ) { dolog=false; };
    mtxconst_fromfile tmpmtx2(col=1)=irnd;
    mtxconst_op tmpmtx2 -= {postexp};
    mtxconst_op tmpmtx2(j) = j^2;
    mtxconst_op postsd += {tmpmtx2};
  };
  mtxconst_op postsd /= N_smp;
  mtxconst_op postsd(j) = sqrt(j);
  calc mtxsum(postsd)/N_el;
  mtxconst_new postsd2(N_el,1,0);
  sfor (i; N_smp) {
    ifstream irnd ( "output/test_25_tmp2_" & {i}$int ) { dolog=false; };
    mtxconst_fromfile tmpmtx2(col=1)=irnd;
    mtxconst_op tmpmtx2 -= {postexp2};
    mtxconst_op tmpmtx2(j) = j^2;
    mtxconst_op postsd2 += {tmpmtx2};
  };
  mtxconst_op postsd2 /= N_smp;
  mtxconst_op postsd2(j) = sqrt(j);
  calc mtxsum(postsd2)/N_el;

sfor (i; N_smp) {
  run_external "rm output/test_25_tmp_" & {i}$int { dolog=false; };
  run_external "rm output/test_25_tmp2_" & {i}$int { dolog=false; };
};

const c = 0;
filter(qid;tmpmtx!{mtxconst_rf_elidvec tmpvec(1,face);}) {
  funplot qid, mtxcoeff(postexp,c,0), mtxcoeff(postexp2,c,0);
  const c = c + 1;
};