123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204 |
- #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;
- }
- 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;
- cout<<"Done"<<endl;
- }
- 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){
- cout<<"Right found"<<endl;
- it_right=it;
- }
- if(b->xright==s.ix){
- cout<<"Left found"<<endl;
- 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->xbottom=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 and hsoil[ix+1]>hsoil[ix]) ++ix;
- if(ix==nX-1)b->position=BorderRight;
- b->xright=ix;
- }
- b->hleft=hsoil[b->xleft];
- b->hright=hsoil[b->xright];
- b->compute_leak_direction();
- // Set hmin
- b->hmin=hsoil[xb];
- basins_to_merge.push_back(b);
- }
- //Determine meta basins`
- 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);
- 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;
- cout<<ix<<" "<<dh<<" "<<runoff_directions[ix]<<endl;
- }
- }
|