~wbetz/fesslix/ssfem

/* Fesslix - Stochastic Analysis
 * Copyright (C) 2010-2017 Wolfgang Betz
 *
 * Fesslix is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Fesslix is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Fesslix.  If not, see <http://www.gnu.org/licenses/>. 
 */


#define VERSION_REV 7
#define VERSION_DATE "2017-05-18"


#include "ssfem_obj.h"
#include "flxlibs_interface.h"


#ifdef __WINDOWS__
  #pragma comment(lib, "flxkernel")
#endif

flxlibs_interface flx_iface;


extern "C" FLXLIBSSFEM_EXPORT void flxlib_init( FlxObjectReadBox* objReadBox, FlxData* dataBox ) {
  GlobalVar.slog(4) << "  SSFEM-library - Version " << VERSION_REV << " (" << VERSION_DATE << ")" << std::endl;
  flx_iface.set_boxes(objReadBox,&(dataBox->FunBox));
  
  flx_iface.add_obj("ssfem", new FlxObjReadSFEM());
  flx_iface.add_obj("ssfemtest", new FlxObjReadSFEMtest()); 
  flx_iface.add_obj("isorfmat", new FlxObjReadIsoRFmat());
  flx_iface.add_obj("isorvmat", new FlxObjReadIsoRVmat());
    
  flx_iface.add_fun("ssfem_moment",new FunReadSFEM_moment() );
  flx_iface.add_fun("ssfem_covariance",new FunReadSFEM_covariance() );
  
  flxSFEM::rbrv_box = &(dataBox->rbrv_box);

}

extern "C" FLXLIBSSFEM_EXPORT void flxlib_close() {
  
  // make sure to unload the SSFEM object
    flxSFEM* sfem = dynamic_cast<flxSFEM*> (FEMbox_get());
    if (sfem) {
      delete sfem;
      FEMbox_get() = NULL;
    }

  flx_iface.free_memory();
}



void FlxObjSFEM_new::task()
{
  if (FEMbox_get()==NULL) {
    throw FlxException("FlxObjSFEM_new::task", "Please load the FEM-library first: 'loadlib \"FEM\";");
  };
  delete FEMbox_get();
  FEMbox_get() = new flxSFEM(p,&data->ConstantBox);
  GlobalVar.slog(4) << "SFEM: created new SFEM calculation" << std::endl;
}

FlxObjReadSFEM::FlxObjReadSFEM() 
{
  // order
    AllDefParaBox->insert(new FlxOptionalParaFun(3,"flxrandom::order"));
    ParaBox.insert("order", "flxrandom::order");
}

int FlxObjReadSFEM::get_order()
{
  FlxOptionalParaBase* p1 = ParaBox.get("order"); 
  void *v1 = p1->get();
  FlxFunction* fun = static_cast<FlxFunction*>( v1 ); 
  int iv = fun->cast2int();
  delete fun;
  if (iv < 0) {
    std::ostringstream ssV_2;
    ssV_2 << "Order of the homogeneous chaos (" << iv << ") must not be a negativ value.";
    throw FlxException("FlxObjReadSFEM::get_order", ssV_2.str(), reader->getCurrentPos() );
  }
  return iv;
}

FlxObjBase* FlxObjReadSFEM::read_new(bool errSerious)
{
  read_optionalPara(errSerious);
  return new FlxObjSFEM_new(get_doLog(),get_order(),get_stream(),get_verbose());
}

FlxObjBase* FlxObjReadSFEMtest::read() {
  //FlxFunction* f1 = new FlxFunction(funReader);
  int M = reader->get_UInt<int>();
  int p = reader->get_UInt<int>();
  read_optionalPara(false);
  return new FlxObjSFEMtest(get_doLog(), M, p, get_stream() );
}


void FlxObjSFEMtest::task()
{
  sout() << "Polynomial Chaos" << std::endl;
  
  PolyChaos* pb = new PolyChaos(M,p);

  tulong N = pb->get_SeqNumb();
  sout() << std::endl;
  sout() << N << ": " << pb->get_ExpPB2(N-1) << std::endl;
  
  delete pb;
}

const std::string FunSFEM_moment::write()
{
  return "sfem_moment(" + write_dof() + "," + moment->write() + ")";
}

const bool FunSFEM_moment::optimize(FunBasePtr& optf, const Fun_OptimizeInfo &foi)
{
  child_optimize(moment,foi);
  return FunFEM_base::optimize(optf,foi);
}

const bool FunSFEM_moment::dependOn_Const(const tdouble*const thenumber)
{
  throw FlxException_NotImplemented("FunSFEM_moment::dependOn_Const");
}

const tdouble FunSFEM_moment::calc()
{
  flxSFEM* sfem = dynamic_cast<flxSFEM*> (FEMbox_get());
  if (sfem == NULL) {
    std::ostringstream ssV;
    ssV << "This is not a stochastic finite element analysis.";
    throw FlxException("FunSFEM_moment::calc_1", ssV.str() );
  }
  
  int nth=int(round_flx(fabs(moment->calc())));
  if (nth < 0) {
    std::ostringstream ssV;
    ssV << "Moment to calculate must not be negative.";
    throw FlxException("FunSFEM_moment::calc_1", ssV.str() );
  } else if ( nth == 0) {
    return ONE;
  }
  
  flxDOFn* dofP = get_dof();
  FlxPeval* sp = sfem->get_Pol(dofP->get_DOFindex());
  tdouble mu = sp->calc_Expectation();
  
  if (nth == 1) {
    delete sp;
    return mu;
  }
  
  mu*=-1;
  sp->add_term(FlxPeval_term(sfem->get_RndBox().get_NRV()),mu);
  
  FlxPeval* spU = sp;
  sp = new FlxPeval(*sp);
  FlxPeval* spT1;
  for ( int i = 2; i <= nth; ++i) {
    spT1 = sp;
    sp = new FlxPeval(*spT1,*spU);
    delete spT1;
  }
  delete spU;
  
  const tdouble cv = sp->calc_Expectation();
  delete sp;
  return cv;
}

FunBase* FunReadSFEM_moment::read(bool errSerious)
{
  FunBase* node;
  FlxDof::dof dofV;
  read_node(node,dofV,errSerious);
  reader->getChar(',',errSerious); 
  FunBase* moment = FunctionList->read(errSerious);
  return new FunSFEM_moment(node,dofV,moment);
}

const std::string FunSFEM_covariance::write()
{
  return "sfem_covariance(" + write_dof() + "," + write_dof(node2,dof2) + ")";
}

const bool FunSFEM_covariance::optimize(FunBasePtr& optf, const Fun_OptimizeInfo& foi)
{
  child_optimize(node2,foi);
  child_optimize(joint2,foi);
  return FunFEM_base::optimize(optf,foi);
}

const bool FunSFEM_covariance::dependOn_Const(const tdouble*const thenumber)
{
  throw FlxException_NotImplemented("FunSFEM_covariance::dependOn_Const");
}

const tdouble FunSFEM_covariance::calc()
{
  flxSFEM* sfem = dynamic_cast<flxSFEM*> (FEMbox_get());
  if (sfem == NULL) {
    std::ostringstream ssV;
    ssV << "This is not a stochastic finite element analysis.";
    throw FlxException("FunSFEM_covariance::calc_1", ssV.str() );
  }
  
  flxDOFn* dofP = get_dof();
  flxDOFn* dofP2 = get_dof(node2,dof2);
  FlxPeval* sp = sfem->get_Pol(dofP->get_DOFindex());
  tdouble mu = sp->calc_Expectation();
  FlxPeval* sp2 = sfem->get_Pol(dofP2->get_DOFindex());
  tdouble mu2 = sp2->calc_Expectation();

  mu*=mu2;
  
  FlxPeval* spT = sp;
  sp = new FlxPeval(*spT,*sp2);
  delete spT; delete sp2;
  mu2 = sp->calc_Expectation();
  delete sp;
  return mu2 - mu;
}

FunBase* FunReadSFEM_covariance::read(bool errSerious)
{
  FunBase* node; FunBase* node2;
  FlxDof::dof dofV; FlxDof::dof dof2;
  read_node(node,dofV,errSerious);
  reader->getChar(',',errSerious);
  read_node(node2,dof2,errSerious);
  return new FunSFEM_covariance(node,dofV,node2,dof2);
}


// -------------------------- IsoRVMaterial ----------------------------------------------------------------------

flxIsoRVMat::flxIsoRVMat(FlxFunction* E, FlxFunction* sD, FlxFunction* nu, FlxFunction* density)
: flxIsotropicMaterial('E',E,'n',nu,density,0),sD(sD)
{

}


flxIsoRVMat::flxIsoRVMat(const flxIsoRVMat& FlxM, const tuint mid)
: flxIsotropicMaterial(FlxM,mid), sD(new FlxFunction(*(FlxM.sD)))
{
  std::ostringstream ssV;
  ssV << "mat_" << mid;
  ts = new RBRV_set_noise_RF(false,1,ssV.str(),false);
  try {
    data->rbrv_box.register_set(ts);
  } catch (FlxException &e) {
    delete ts;
    ts = NULL;
    throw;
  }
  y_smp = &(ts->get_y()[0]);
}

flxIsoRVMat* flxIsoRVMat::copy(const tuint mid)
{
  return new flxIsoRVMat(*this,mid);
}

void flxIsoRVMat::log_info(std::ostream& lout)
{
  lout << "  isorvmat: ";
  lout << v1c << "=" << v1->write() << "; " << v2c << "=" << v2->write() << "; ";
  lout << "Random in E: ";
  lout << "Normal(mu=" << v1->write() << "; sd=" << sD->write() << "); ";
  lout << std::endl;
}

const tdouble flxIsoRVMat::get_E() const
{
  const int Current_rv_index = current_RF_mode_get();
  if (Current_rv_index < 0) {
    tdouble d = flxIsotropicMaterial::get_E();
    d+=sD->cast2positive()*(*y_smp);
    return d;
  } else if (Current_rv_index == 0) {
    return flxIsotropicMaterial::get_E();
  } else if (static_cast<tuint>(Current_rv_index)-1 == ts->get_RndBox_index()) {
    return sD->cast2positive();
  } else return ZERO;
}

const tdouble flxIsoRVMat::get_G() const
{
  const tdouble E = get_E();
  const tdouble nu = get_nu();
  return check_G(E/(2*(ONE+nu)));
}

const tdouble flxIsoRVMat::get_lambda() const
{
  const tdouble E = get_E();
  const tdouble nu = get_nu();
  return E*nu/((ONE+nu)*(ONE-2*nu));
}

FlxFunction* FlxObjReadIsoRVmat::read_matprop(char& c, bool errSerious)
{
  c = 0;
  const std::string s1 = reader->getWord(true,errSerious);
  if (s1 == "e") {
    c = 'E';
  } else if (s1 == "sd") {
    c = 's';
  } else if (s1 == "nu") {
    c = 'n';
  } else { // ERROR: not known !!!
    std::ostringstream ssV_2;
    ssV_2 << "Material parameter '" << s1 << "' is not known.";
    FlxError(errSerious,"FlxObjReadIsoRVmat::read_matprop_1", ssV_2.str(), reader->getCurrentPos() );
  }
  return new FlxFunction(funReader,errSerious);
}

void FlxObjReadIsoRVmat::get_matprop(FlxFunctionPtr& E, FlxFunctionPtr& sD, FlxFunctionPtr& nu, bool errSerious)
{
  bool b1=true;
  do {
    char c1=0;
    FlxFunction* f=read_matprop(c1,errSerious);
    if (c1=='E') {
      if (E!=NULL) {
	std::ostringstream ssV_2;
	ssV_2 << "'E' was already defined.";
	FlxError(errSerious,"FlxObjReadIsoRVmat::read_1", ssV_2.str(), reader->getCurrentPos() );
      }
      E=f;
    } else if (c1=='s') {
      if (sD!=NULL) {
	std::ostringstream ssV_2;
	ssV_2 << "'sD' was already defined.";
	FlxError(errSerious,"FlxObjReadIsoRVmat::read_2", ssV_2.str(), reader->getCurrentPos() );
      }
      sD=f;
    } else if (c1=='n') {
      if (nu!=NULL) {
	std::ostringstream ssV_2;
	ssV_2 << "'nu' was already defined.";
	FlxError(errSerious,"FlxObjReadIsoRVmat::read_4", ssV_2.str(), reader->getCurrentPos() );
      }
      nu=f;
    }
    if (E!=NULL && sD!=NULL && nu!=NULL) b1=false;
  } while (b1);
}

FlxObjBase* FlxObjReadIsoRVmat::read()
{
  FlxFunction* index = new FlxFunction(funReader,false);
  FlxFunction *E=NULL, *sD=NULL, *nu=NULL;
  try {
    get_matprop(E,sD,nu);
    read_optionalPara(false);
    flxIsoRVMat* mat = new flxIsoRVMat(E,sD,nu,get_optPara_FlxFunction("density"));
    FlxObjBase* obj1 = new FlxObjIsoMaterial(get_doLog(),index, mat);
    return obj1;  
  } catch (FlxException &e) {
    FLXMSG("FlxObjReadIsoRVmat::read",1);
    delete index;
    if(E) delete E;
    if(sD) delete sD;
    if(nu) delete nu;
    throw;
  }
}

// -------------------------- IsoRFMaterial ----------------------------------------------------------------------


flxIsoRFMat::flxIsoRFMat(FlxRandomField_base* RFptr, FlxFunction* nu, FlxFunction* density, const tuint deg_E)
: flxIsotropicMaterial('E',RFptr->get_mean_fun(),'n',nu,density,deg_E), RFptr(RFptr)
{
  FlxRandomField_KL_Analytical_ExpCorr_Gauss_1D* klanalytgauss1d = dynamic_cast<FlxRandomField_KL_Analytical_ExpCorr_Gauss_1D*>(RFptr);
  if (klanalytgauss1d==NULL) {
    FlxRandomField_KL_FEM_Gauss* klfemgauss = dynamic_cast<FlxRandomField_KL_FEM_Gauss*>(RFptr);
    if (klfemgauss==NULL) {
      std::ostringstream ssV;
      ssV << "The specified random field is not a Gaussian random field.";
      throw FlxException("flxIsoRFMat::flxIsoRFMat", ssV.str() );
    }
  }
}

flxIsoRFMat::flxIsoRFMat(const flxIsoRFMat& FlxM, const tuint mid): flxIsotropicMaterial(FlxM,mid), RFptr(FlxM.RFptr)
{
  RFptr->set_automatic();
}

flxIsoRFMat* flxIsoRFMat::copy(const tuint mid)
{
  return new flxIsoRFMat(*this,mid);
}

void flxIsoRFMat::log_info(std::ostream& lout)
{
  lout << "  isorfmat: ";
  lout << v1c << "=" << v1->write() << "; " << v2c << "=" << v2->write() << "; ";
  lout << "Random in E: ";
  RFptr->print_para(lout);
  lout << std::endl;
}

const bool flxIsoRFMat::depend_E_on_Const(const tdouble* const thenumber) {
  if (ConstantBox->is_SpatialVar(thenumber)) return true;
  else return flxIsotropicMaterial::depend_E_on_Const(thenumber);
}

const tdouble flxIsoRFMat::get_E() const
{
  const int Current_rv_index = current_RF_mode_get();
  if (!RFptr->is_Assembled() || Current_rv_index==0) return flxIsotropicMaterial::get_E();
  if (Current_rv_index == -1) {
    return RFptr->get_realization();
  } else {
    return RFptr->calc_Hi(Current_rv_index,true);
  }
}

const tdouble flxIsoRFMat::get_G() const
{
  const tdouble E = get_E();
  const tdouble nu = get_nu();
  return check_G(E/(2*(ONE+nu)));
}

const tdouble flxIsoRFMat::get_lambda() const
{
  const tdouble E = get_E();
  const tdouble nu = get_nu();
  return E*nu/((ONE+nu)*(ONE-2*nu));
}

FlxObjBase* FlxObjReadIsoRFmat::read()
{
  FlxFunction* index = new FlxFunction(funReader);
  FlxFunction *RFindex=NULL,*nu=NULL;
  try {
    RFindex = new FlxFunction(funReader);
    reader->getWord("nu",false);
    nu = new FlxFunction(funReader);
    read_optionalPara(false);
    FlxRandomField_base* RFptr = RFbox_get().get(RFindex->cast2tuintW0(false)); 
    delete RFindex; RFindex = NULL;
    RFindex = get_optPara_FlxFunction("degreee");
    const tuint deg_E = RFindex->cast2tuintW0(false);
    delete RFindex; RFindex = NULL;
    flxIsoRFMat *mat = new flxIsoRFMat(RFptr,nu,get_optPara_FlxFunction("density"),deg_E);
    return new FlxObjIsoMaterial(get_doLog(), index, mat);
  } catch ( FlxException& e) {
    FLXMSG("FlxObjReadIsoRFmat::read",1);
    delete index;
    if (RFindex) delete RFindex;
    if (nu) delete nu;
    throw;
  }
}