#include "set_layers.h" #define ind(r,c) ((r)*N + (c)) void rand1d(double * arr, int N) { if( N == 1 ) return; const double t = 0.15; double a = 0.5+(2*(rand()*1.0/RAND_MAX)-1)*t; // in 0.5+[-1,1]*t arr[N/2] = arr[0]*a + arr[N]*(1-a); rand1d(arr,N/2); rand1d(arr+N/2,N-N/2); } std::pair layer1d(double * arr, int N, double z) { for(int k = 0; k < N; ++k) { if( z >= arr[k] && z <= arr[k+1] ) return std::make_pair(k, (z-arr[k])/(arr[k+1]-arr[k])); } return std::make_pair(-1,0.0); } const double noind = std::numeric_limits::quiet_NaN(); void rand2d(double * arr, int N, int Nl, int Nr, int Nb, int Nt, double t) { //std::cout << "Nl:Nr " << Nl <<":" << Nr << " Nb:Nt " << Nb << ":" << Nt << std::endl; if( Nr - Nl < 2 && Nt - Nb < 2 ) { //std::cout << "exit" << std::endl; return; } //const double t = 0.15; int Nk = (Nb+Nt)/2; int Nm = (Nl+Nr)/2; //std::cout << "Nk " << Nk << " Nm " << Nm << std::endl; double lb = arr[ind(Nb,Nl)]; double rb = arr[ind(Nb,Nr)]; double lt = arr[ind(Nt,Nl)]; double rt = arr[ind(Nt,Nr)]; if( lb != lb || rb != rb || lt != lt || rt != rt ) throw -1; if( arr[ind(Nk,Nl)] != arr[ind(Nk,Nl)] ) arr[ind(Nk,Nl)] = 0.5*(lb + lt) + (2*(rand()*1.0/RAND_MAX)-1)*t; if( arr[ind(Nk,Nr)] != arr[ind(Nk,Nr)] ) arr[ind(Nk,Nr)] = 0.5*(rb + rt) + (2*(rand()*1.0/RAND_MAX)-1)*t; if( arr[ind(Nb,Nm)] != arr[ind(Nb,Nm)] ) arr[ind(Nb,Nm)] = 0.5*(lb + rb) + (2*(rand()*1.0/RAND_MAX)-1)*t; if( arr[ind(Nt,Nm)] != arr[ind(Nt,Nm)] ) arr[ind(Nt,Nm)] = 0.5*(lt + rt) + (2*(rand()*1.0/RAND_MAX)-1)*t; arr[ind(Nk,Nm)] = 0.25*(lb+rb+lt+rt) + (2*(rand()*1.0/RAND_MAX)-1)*t; rand2d(arr,N,Nl,Nm,Nb,Nk,t*0.5); rand2d(arr,N,Nm,Nr,Nb,Nk,t*0.5); rand2d(arr,N,Nl,Nm,Nk,Nt,t*0.5); rand2d(arr,N,Nm,Nr,Nk,Nt,t*0.5); } void init2d(double * arr, int N, double mint, double maxt) { for(int k = 0; k < N*N; ++k) arr[k] = noind; arr[ind(0 ,0 )] = mint+(rand()*1.0/RAND_MAX)*(maxt-mint); arr[ind(0 ,N-1)] = mint+(rand()*1.0/RAND_MAX)*(maxt-mint); arr[ind(N-1,0 )] = mint+(rand()*1.0/RAND_MAX)*(maxt-mint); arr[ind(N-1,N-1)] = mint+(rand()*1.0/RAND_MAX)*(maxt-mint); } double intrp2d(double * arr, int N, double x, double y) { int n = ceil(x*(N-1)); int m = ceil(y*(N-1)); if( n == 0 ) n = 1; if( m == 0 ) m = 1; double dh = 1.0/(double)(N-1); double kx = (x-(n-1)*dh)/dh; double ky = (y-(m-1)*dh)/dh; //if( kx < 0 || kx > 1 ) std::cout << "bad kx: " << kx << " x is " << x << " n is " << n << " dh is " << dh << " N is " << N << std::endl; //if( ky < 0 || ky > 1 ) std::cout << "bad ky: " << ky << " y is " << y << " m is " << m << " dh is " << dh << " N is " << N << std::endl; double lb = arr[ind(m-1,n-1)]; double rb = arr[ind(m-1,n+0)]; double lt = arr[ind(m+0,n-1)]; double rt = arr[ind(m+0,n+0)]; if( lb != lb || rb != rb || lt != lt || rt != rt ) throw -1; return (1-ky)*(lb*(1-kx) + rb*kx) + ky*(lt*(1-kx) + rt*kx); } void SetLayers::RandomLayersZ(int num_layers) { layers_z.resize(num_layers+1); layers_z[0] = 0; layers_z[num_layers] = 1; rand1d(&layers_z[0],num_layers); } void SetLayers::DeformHeight(double height) { const int N = 128; std::vector map(N*N,0.0); init2d(&map[0],N,0.0,height); rand2d(&map[0],N,0,N-1,0,N-1,height*0.5); for(int k = 0; k < N*N; ++k) if( map[k] != map[k] ) std::cout << "NaN at row " << k/N << " col " << k%N << std::endl; double cmin[3], cmax[3]; for(int d = 0; d < 3; ++d) { cmin[d]=1.0e+20; cmax[d]=-1.0e+20; } for(Mesh::iteratorNode n = mesh.BeginNode(); n != mesh.EndNode(); ++n) { for(int d = 0; d < 3; ++d) { if( cmin[d] > n->Coords()[d] ) cmin[d] = n->Coords()[d]; if( cmax[d] < n->Coords()[d] ) cmax[d] = n->Coords()[d]; } } for(Mesh::iteratorNode n = mesh.BeginNode(); n != mesh.EndNode(); ++n) { double c[3] = {0,0,0}; for(int d = 0; d < 3; ++d) c[d] = (n->Coords()[d]-cmin[d])/(cmax[d]-cmin[d]); n->Coords()[2] += intrp2d(&map[0],N,c[0],c[1]); } } void SetLayers::DeformLayers(double coef) { const int N = 128; std::vector< std::vector > map_layer(layers_z.size()); #if defined(USE_OMP) #pragma omp parallel for #endif for(int k = 0; k < map_layer.size(); ++k) { map_layer[k].resize(N*N); init2d(&map_layer[k][0],N,1-coef,1); rand2d(&map_layer[k][0],N,0,N-1,0,N-1,coef); for(int kk = 0; kk < N*N; ++kk) if( map_layer[k][kk] < 0.0 ) map_layer[k][kk] = 0; } double cmin[3], cmax[3]; for(int d = 0; d < 3; ++d) { cmin[d]=1.0e+20; cmax[d]=-1.0e+20; } for(Mesh::iteratorNode n = mesh.BeginNode(); n != mesh.EndNode(); ++n) { for(int d = 0; d < 3; ++d) { if( cmin[d] > n->Coords()[d] ) cmin[d] = n->Coords()[d]; if( cmax[d] < n->Coords()[d] ) cmax[d] = n->Coords()[d]; } } TagInteger layer_tag = mesh.CreateTag("LAYER",DATA_INTEGER,CELL,NONE,1); TagReal coef_tag = mesh.CreateTag("LAYER_COEF",DATA_REAL,CELL,NONE,1); for(Mesh::iteratorCell c = mesh.BeginCell(); c != mesh.EndCell(); ++c) { double cnt[3], cc[3] = {0,0,0}; c->Centroid(cnt); for(int d = 0; d < 3; ++d) cc[d] = (cnt[d]-cmin[d])/(cmax[d]-cmin[d]); std::pair lc = layer1d(&layers_z[0],layers_z.size()-1,cc[2]); layer_tag[*c] = lc.first; coef_tag[*c] = lc.second; } for(Mesh::iteratorNode n = mesh.BeginNode(); n != mesh.EndNode(); ++n) { double c[3] = {0,0,0}; for(int d = 0; d < 3; ++d) c[d] = (n->Coords()[d]-cmin[d])/(cmax[d]-cmin[d]); std::pair lc = layer1d(&layers_z[0],layers_z.size()-1,c[2]); if( lc.first == -1 ) std::cout << "layer not found for " << c[0] << "," << c[1] << "," << c[2] << " node " << n->LocalID() << std::endl; double h = 0; int kk = lc.first; for(int k = 0; k < kk; ++k) h += (layers_z[k+1]-layers_z[k])*(cmax[2]-cmin[2])*intrp2d(&map_layer[k][0],N,c[0],c[1]); h += lc.second*(layers_z[kk+1]-layers_z[kk])*(cmax[2]-cmin[2])*intrp2d(&map_layer[kk][0],N,c[0],c[1]); n->Coords()[2] = h; } //TODO: fix zero edges int bad_edge = 0,collapse_edge =0; mesh.BeginModification(); for(Mesh::iteratorEdge e = mesh.BeginEdge(); e != mesh.EndEdge(); ++e) if( !e->Hidden() ) { if( e->Length() < 1.0e-9 ) { if( !e->Hidden() ) { Edge::CollapseEdge(e->self(),0); collapse_edge++; } bad_edge++; } } mesh.EndModification(); std::cout << "Bad edges: " << bad_edge << " collapsed " << collapse_edge << std::endl; }