RichardsFastSlow/src/kernel/geometry.cpp

#include "geometry.hpp"
#include "physics.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;
}


Geometry::~Geometry(){
  if(hsoil!=nullptr){
    delete[] hsoil;
    delete[] dhsoil;
    delete[] hbot;
    delete[] dhbot;
    delete[] dZ;
    delete[] nZ;
    for(size_t i=0;i<nX;++i){
      delete[] Z[i];
    }
    delete[] Z;
  }
}

void
Geometry::initZ(double dZ_avg,bool init){
  if(init){
    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+1;
    if(init) Z[k]=new double[n+1];
    for(size_t j=0;j<n+1;++j){
      Z[k][j]=hbot[k]+j*dz;
    }
  }
}

void
Geometry::save(fstream& file) {
  file.write((char*)&lX,sizeof(double));
  file.write((char*)&nX,sizeof(size_t));
  file.write((char*)&dX,sizeof(double));
  file.write((char*)hsoil,nX*sizeof(double));
  file.write((char*)dhsoil,nX*sizeof(double));
  file.write((char*)hbot,nX*sizeof(double));
  file.write((char*)dhbot,nX*sizeof(double));
  file.write((char*)nZ,nX*sizeof(size_t));
  file.write((char*)dZ,nX*sizeof(double));
  for(size_t i=0;i<nX;++i){
    file.write((char*)Z[i],nZ[i]*sizeof(double));
  }
}

void
Geometry::load(fstream& file,bool init) {
  file.read((char*)&lX,sizeof(double));
  file.read((char*)&nX,sizeof(size_t));
  file.read((char*)&dX,sizeof(double));
  if(init){
    hsoil=new double[nX];
    dhsoil=new double[nX];
    hbot=new double[nX];
    dhbot=new double[nX];
    nZ=new size_t[nX];
    dZ=new double[nX];
    Z=new double*[nX];
  }
  file.read((char*)hsoil,nX*sizeof(double));
  file.read((char*)dhsoil,nX*sizeof(double));
  file.read((char*)hbot,nX*sizeof(double));
  file.read((char*)dhbot,nX*sizeof(double));
  file.read((char*)nZ,nX*sizeof(size_t));
  file.read((char*)dZ,nX*sizeof(double));
  for(size_t i=0;i<nX;++i){
    if(init) Z[i]=new double[nZ[i]];
    file.read((char*)Z[i],nZ[i]*sizeof(double));
  }
  //Comoute dhsoil and dhbot
  dhsoil[0]=(hsoil[1]-hsoil[0])/dX;
  dhbot[0]=(hbot[1]-hbot[0])/dX;

  for(size_t i=1;i<nX-1;++i){
    dhsoil[i]=(hsoil[i+1]-hsoil[i-1])/(2*dX);
    dhbot[i]=(hbot[i+1]-hbot[i-1])/(2*dX);
  }
  dhsoil[nX-1]=(hsoil[nX-1]-hsoil[nX-2])/dX;
  dhbot[nX-1]=(hbot[nX-1]-hbot[nX-2])/dX;

}

pair<list<Basin*>::iterator,list<Basin*>::iterator>
Geometry::find_basins(const Summit& s,list<Basin*>& basins){
  list<Basin*>::iterator it_left;
  list<Basin*>::iterator it_right;
  for(list<Basin*>::iterator it=basins.begin();it!=basins.end();++it){
    Basin* b=*it;
    if(b->xleft==s.ix){
      it_right=it;
    }
    if(b->xright==s.ix){
      it_left=it;
    }
  }
  return pair<list<Basin*>::iterator,list<Basin*>::iterator>(it_left,it_right);
}

void Geometry::compute_basins(){
  //Find local min point of hsoil
  vector<size_t> v_min;
  if(hsoil[0]<hsoil[1]) v_min.push_back(0);
  for(size_t ix=1;ix<nX-1;++ix){
    double h=hsoil[ix];
    if(hsoil[ix-1]>h and h<hsoil[ix+1]) v_min.push_back(ix);
  }
  if(hsoil[nX-2]>hsoil[nX-1]) v_min.push_back(nX-1);
  size_t nb=v_min.size();

  //Primitive Basins
  vector<Summit> summits;
  list<Basin*> basins_to_merge;

  for(size_t ib=0;ib<nb;++ib){
    Basin* b=new Basin;
    size_t xb=v_min[ib];
    //Find left bound
    b->xacc=xb;
    if(xb==0){
      b->xleft=xb;
      b->position=BorderLeft;
    }
    else{
      size_t ix=xb;
      while(ix>0 and hsoil[ix-1]>hsoil[ix]) --ix;
      b->xleft=ix;
      if(ix==0) b->position=BorderLeft;
      if(ix>0) summits.emplace_back(ix,hsoil[ix]);
    }
    //Find right bound
    if(xb==nX-1){
      b->xright=xb;
      b->position=BorderRight;
    }
    else{
      size_t ix=xb;
      while(ix<nX-1 and hsoil[ix+1]>hsoil[ix]) ++ix;
      if(ix==nX-1)b->position=BorderRight;
      b->xright=ix;
    }
    b->hleft=hsoil[b->xleft]+Physics::minimal_height_to_runoff;
    b->hright=hsoil[b->xright]+Physics::minimal_height_to_runoff;

    b->compute_leak_direction();
    b->compute_maxs(hsoil);

    // Set hmin
    b->hmin=hsoil[xb]+Physics::minimal_height_to_runoff;

    basins_to_merge.push_back(b);
  }
  //Determine meta basins

  if(summits.empty()){
    root_basin=new Basin;
    root_basin->position=BorderBoth;
    root_basin->left=nullptr;
    root_basin->right=nullptr;
    root_basin->xacc=0;
    root_basin->leak_direction=None;
    root_basin->xleft=0;
    root_basin->xright=nX-1;
    double h=hsoil[0];
    root_basin->hleft=h;
    root_basin->hright=h;
    root_basin->hmin=h;

  }
  else{
    sort(summits.begin(),summits.end());
    for(auto it=summits.begin();it!=summits.end();++it){
      Summit& s=*it;
      auto res=find_basins(s,basins_to_merge);
      list<Basin*>::iterator it_left=res.first;
      list<Basin*>::iterator it_right=res.second;
      Basin* b_left=*it_left;
      Basin* b_right=*it_right;

      basins_to_merge.erase(it_left);
      basins_to_merge.erase(it_right);
      Basin* b=new Basin(b_left,b_right);
      b->compute_maxs(hsoil);
      basins_to_merge.push_back(b);
    }
    root_basin=basins_to_merge.front();
    root_basin->display();
  }

  // Compute Runoff directions
  runoff_directions=new Direction[nX];
  for(size_t ix=1;ix<nX-1;++ix){
    double dh=(hsoil[ix+1]-hsoil[ix-1])/(2*dX);
    if(abs(dh)<max_slope_both_side_runoff) runoff_directions[ix]=Both;
    else runoff_directions[ix]=(dh>0)?Left:Right;
  }
  for(size_t ib=0;ib<nb;++ib){
    size_t xb=v_min[ib];
    runoff_directions[xb]=None;
  }
  runoff_directions[0]=Right;
  runoff_directions[nX-1]=Left;
}