RichardsFastSlow/src/geometry.cpp

#include "geometry.hpp"

double inf = numeric_limits<double>::infinity();

Geometry::Geometry(){
  lX=0;
  nX=0;
  dX=0;
  hsoil=nullptr;
  dhsoil=nullptr;
  hbot=nullptr;
  dhbot=nullptr;
  nZ=nullptr;
  dZ=nullptr;
  Z=nullptr;
}

// Obsolete, we must add hsat
Geometry::Geometry(double _lX,size_t _nX,size_t _nZ,Func _hsoil,Func _dhsoil,Func _hbot,Func _dhbot){
  lX=_lX;
  nX=_nX;
  dX=lX/nX;
  hsoil=new double[nX];
  dhsoil=new double[nX];
  hbot=new double[nX];
  dhbot=new double[nX];
  double hs_max=-inf,hb_min=inf;
  double v;
  for(size_t k=0;k<nX;++k){
    double x=k*dX;
    v=hsoil[k]=_hsoil(x);
    hs_max=max(v,hs_max);
    dhsoil[k]=_dhsoil(x);
    v=hbot[k]=_hbot(x);
    hb_min=min(v,hb_min);
    dhbot[k]=_dhbot(x);
  }
  double dZ_avg=(hs_max-hb_min)/_nZ;
  initZ(dZ_avg);
}


/*Geometry::Geometry(double _lX,size_t _nX,double depth,size_t _nZ,Spline& Ssoil,Spline& Sbot){
  lX=_lX;
  nX=_nX;
  dX=lX/nX;
  hsoil=new double[nX];
  dhsoil=new double[nX];
  hbot=new double[nX];
  dhbot=new double[nX];
  double hs_max=-inf,hb_min=inf;
  double v;
  for(size_t k=0;k<nX;++k){
    double x=k*dX/lX;
    v=hsoil[k]=1-Ssoil(x);
    hs_max=max(v,hs_max);
    dhsoil[k]=-Ssoil.derivate(x);
    v=hbot[k]=1-Sbot(x);
    hb_min=min(v,hb_min);
    dhbot[k]=-Sbot.derivate(x);
  }
  double dZ_avg=(hs_max-hb_min)/_nZ;
  initZ(dZ_avg);
  }*/

void
Geometry::initZ(double dZ_avg){
  nZ=new size_t[nX];
  dZ=new double[nX];
  Z=new double*[nX];
  for(size_t k=0;k<nX;++k){
    double d=hsoil[k]-hbot[k];
    size_t n=d/dZ_avg;
    double dz=d/n;
    dZ[k]=dz;
    nZ[k]=n;
    Z[k]=new double[n];
    for(size_t j=0;j<n;++j){
      Z[k][j]=hbot[k]+j*dz;
    }
  }
}

void
Geometry::initDefault(){
  lX=20;
  nX=100;
  dX=lX/(nX-1);
  hsoil=new double[nX];
  dhsoil=new double[nX];
  hbot=new double[nX];
  dhbot=new double[nX];
  for(size_t i=0;i<nX;++i){
    hsoil[i]=0.7;
    hbot[i]=0.3;
    dhsoil[i]=0;
    dhbot[i]=0;
  }
  depth=0.7-0.3;
  nZ_max=100;
  initZ(depth/nZ_max);
}

void
Geometry::update(double _lX,size_t _nX,double _depth,size_t _nZ_max){
  lX=_lX;
  depth=_depth;
  if(_nX!=nX or _nZ_max!=nZ_max){
    delete[] hsoil;
    delete[] dhsoil;
    delete[] hbot;
    delete[] dhbot;
    delete[] nZ;
    delete[] dZ;
    for(size_t i=0;i<nX;++i) delete[] Z[i];
    delete[] Z;
    nX=_nX;
    nZ_max=_nZ_max;
    hsoil=new double[nX];
    dhsoil=new double[nX];
    hbot=new double[nX];
    dhbot=new double[nX];
    initZ(depth/nZ_max);
  }
}