#include "Temparray.h" Temparray::Temparray(float initialtemp, int x0, int y0, int z0){ temparraynew = new float[x0*y0*z0*6*4]; temparrayold = new float[x0*y0*z0*6*4]; cubearray = new Cubehole[x0*y0*z0*6]; weight1 = new int[x0*y0*z0]; sx = x0; sy = y0; sz = z0; tempInit(initialtemp, x0, y0, z0); // static const float pos[5] = {-2.0, -1.0, 0.0, 1.0, 2.0}; for(int i = 0; i < x0; ++i) { for(int j = 0; j < y0; ++j) { for(int k = 0; k < z0; ++k) { for(int l = 0; l < 6; ++l) { cubehole(i, j, k, l).setSize((6-l)/6.0*0.9, 0.15, (6-l)/6.0*0.9, (5-l)/6.0*0.9, (5-l)/6.0*0.9); float x, y, z; if(x0 % 2 == 0) x = -(x0/2)+i+0.5; else if(x0 % 2 == 1) x = -(x0-1)/2+i; if(y0 % 2 == 0) y = -(y0/2)+j+0.5; else if(y0 % 2 == 1) y = -(y0-1)/2+j; if(z0 % 2 == 0) z = -(z0/2)+k+0.5; else if(z0 % 2 == 1) z = -(z0-1)/2+k; cubehole(i, j, k, l).setPos(x, y, z); cubehole(i, j, k, l).setColor(vmml::vec4f(1.0, 0.5, 1.0, 1.0), vmml::vec4f(0.5, 0.5, 1.0, 1.0), vmml::vec4f(1.0, 0.5, 0.0, 1.0), vmml::vec4f(0.0, 0.5, 0.0, 1.0)); } } } } calcweight(); } void Temparray::calcTemp(){ float conductivity = 0.7, specificcapacity = 1513; float areasmall, area, areabig, distance, capacity, capacity2, volume, thermalresistance, thermalresistance2, tau12, tau21; float meterperunit = 3.0; float width, height, depth, innerwidth, innerdepth; float width2, height2, depth2, innerwidth2, innerdepth2; float loadableenergy, unloadableenergy; calcweight(); for(int i = 0; i < sx; ++i) { for(int j = 0; j < sy; ++j) { for(int k = 0; k < sz; ++k) { for(int l = 1; l < 6; ++l) { for(int m = 0; m < 4; ++m) { width = cubehole(i, j, k, l).getWidth(); width2 = cubehole(i, j, k, l-1).getWidth(); height = cubehole(i, j, k, l).getHeight(); height2 = cubehole(i, j, k, l-1).getHeight(); depth = cubehole(i, j, k, l).getDepth(); depth2 = cubehole(i, j, k, l-1).getDepth(); innerwidth = cubehole(i, j, k, l).getInnerWidth(); innerwidth2 = cubehole(i, j, k, l-1).getInnerWidth(); innerdepth = cubehole(i, j, k, l).getInnerDepth(); innerdepth2 = cubehole(i, j, k, l-1).getInnerDepth(); if(m % 2 == 0) { area = (((width + innerwidth)/2)*height) / 0.9 * meterperunit; areabig = (((width2 + innerwidth2)/2)*height2) / 0.9 * meterperunit; thermalresistance = ((depth2 / 2 - (depth2 - innerdepth2)/4) - (depth / 2 - (depth - innerdepth)/4)) / 0.9 * meterperunit / (conductivity * areabig); thermalresistance2 = ((depth2 / 2 - (depth2 - innerdepth2)/4) - (depth / 2 - (depth - innerdepth)/4)) / 0.9 * meterperunit / (conductivity * area); capacity = specificcapacity *((((depth - innerdepth)/2)*((width + innerwidth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((depth2 - innerdepth2)/2)*((width2 + innerwidth2)/2)* height2)/0.9*meterperunit); } else if(m % 2 == 1) { area = (((depth + innerdepth)/2)*height) / 0.9 * meterperunit; areabig = (((depth2 + innerdepth2)/2)*height2) / 0.9 * meterperunit; thermalresistance = ((width2 / 2 - (width2 - innerwidth2)/4) - (width / 2 - (width - innerwidth)/4)) / 0.9 * meterperunit / (conductivity * areabig); thermalresistance2 = ((width2 / 2 - (width2 - innerwidth2)/4) - (width / 2 - (width - innerwidth)/4)) / 0.9 * meterperunit / (conductivity * area); capacity = specificcapacity *((((depth + innerdepth)/2)*((width - innerwidth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((depth2 + innerdepth2)/2)*((width2 - innerwidth2)/2)* height2)/0.9*meterperunit); } tau12 = capacity * thermalresistance; tau21 = capacity2 * thermalresistance2; temperaturenew(i, j, k, l, m) = temperaturenew(i, j, k, l, m) - ((temperatureold(i, j, k, l, m) - temperatureold(i, j, k, l-1, m))*(1-exp(-(1/tau12)))); temperaturenew(i, j, k, l-1, m) = temperaturenew(i, j, k, l-1, m) + ((temperatureold(i, j, k, l, m) - temperatureold(i, j, k, l-1, m))*(1-exp(-(1/tau21)))); } } } } } for(int h = 0; h <= maxweight; ++h) { for(int i = 0; i < sx; ++i) { for(int j = 0; j < sy; ++j) { for(int k = 0; k < sz; ++k) { for(int l = 0; l < 6; ++l) { for(int m = 0; m < 4; m++) { width = cubehole(i, j, k, l).getWidth(); height = cubehole(i, j, k, l).getHeight(); depth = cubehole(i, j, k, l).getDepth(); innerwidth = cubehole(i, j, k, l).getInnerWidth(); innerdepth = cubehole(i, j, k, l).getInnerDepth(); if(m % 2 == 0)area = (((depth - innerdepth)/2)*((width + innerwidth)/2)) / 0.9 * meterperunit; if(m % 2 == 1)area = (((width - innerwidth)/2)*((depth + innerdepth)/2)) / 0.9 * meterperunit; if(m % 2 == 0)thermalresistance = ((depth - innerdepth)/2) /0.9 * meterperunit /(conductivity * area); if(m % 2 == 1)thermalresistance = ((width - innerwidth)/2) /0.9 * meterperunit /(conductivity * area); if(m % 2 == 0)capacity = specificcapacity *((((depth - innerdepth)/2)* ((width + innerwidth)/2)*height) / 0.9 * meterperunit); if(m % 2 == 1)capacity = specificcapacity *((((width - innerwidth)/2)* ((depth + innerdepth)/2)*height) / 0.9 * meterperunit); loadableenergy = (probetemp - temperatureold(i, j, k, l, m))/thermalresistance/10; unloadableenergy = (temperatureold(i, j, k, l, m)-21.5)/thermalresistance/10; tau12 = capacity * thermalresistance; if(input>=0 && weight(i, j, k)==h){ if(loadableenergy <= 0) inputpower = 0; else if (input >= loadableenergy) { inputpower = loadableenergy; input -= loadableenergy; } else if (input < loadableenergy) { inputpower = input; input = 0; } if(inputpower > 0) { factor = inputpower/loadableenergy; temperaturenew(i, j, k, l, m) = temperaturenew(i, j, k, l, m) + (probetemp-temperatureold(i, j, k, l, m)) *thermalresistance*factor*(1-exp(-(1/tau12))); } } if (input < 0 && maxweight-weight(i, j, k)==h){ if(unloadableenergy <= 0) inputpower = 0; else if (input > unloadableenergy) { inputpower = unloadableenergy; input -= unloadableenergy; } else if (input <= unloadableenergy) { inputpower = input; input = 0; } if(inputpower > 0) { factor = inputpower/unloadableenergy; temperaturenew(i, j, k, l, m) = temperaturenew(i, j, k, l, m) + (21.5-temperatureold(i, j, k, l, m)) *thermalresistance*factor*(1-exp(-(1/tau12))); } } } } } } } } for(int i = 0; i < sx; ++i) { for(int j = 0; j < sy; ++j) { for(int k = 0; k < sz; ++k) { for(int l = 0; l < 6; ++l) { for(int m = 0; m < 4; m++) { width = cubehole(i, j, k, l).getWidth(); height = cubehole(i, j, k, l).getHeight(); depth = cubehole(i, j, k, l).getDepth(); innerwidth = cubehole(i, j, k, l).getInnerWidth(); innerdepth = cubehole(i, j, k, l).getInnerDepth(); if(j > 0) { width2 = cubehole(i, j-1, k, l).getWidth(); height2 = cubehole(i, j-1, k, l).getHeight(); depth2 = cubehole(i, j-1, k, l).getDepth(); innerwidth2 = cubehole(i, j-1, k, l).getInnerWidth(); innerdepth2 = cubehole(i, j-1, k, l).getInnerDepth(); if(m % 2 == 0) { area = (((width + innerwidth)/2) * ((depth-innerdepth)/2)) / 0.9 * meterperunit; areabig = (((width2 + innerwidth2)/2) * ((depth2-innerdepth2)/2)) / 0.9 * meterperunit; thermalresistance = (height/2+height2/2) / 0.9 * meterperunit / (conductivity * ((area + areabig) / 2)); capacity = specificcapacity *((((depth - innerdepth)/2)*((width + innerwidth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((depth2 - innerdepth2)/2)*((width2 + innerwidth2)/2)* height2)/0.9*meterperunit); } else if(m % 2 == 1) { area = (((depth + innerdepth)/2) * ((width-innerwidth)/2)) / 0.9 * meterperunit; areabig = (((depth2 + innerdepth2)/2) * ((width2-innerwidth2)/2)) / 0.9 * meterperunit; thermalresistance = (height/2+height2/2) / 0.9 * meterperunit / (conductivity * ((area + areabig) / 2)); capacity = specificcapacity *((((depth + innerdepth)/2)*((width - innerwidth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((depth2 + innerdepth2)/2)*((width2 - innerwidth2)/2)* height2)/0.9*meterperunit); } tau12 = capacity * thermalresistance; tau21 = capacity2 * thermalresistance; temperaturenew(i, j, k, l, m) = temperaturenew(i, j, k, l, m) - ((temperatureold(i, j, k, l, m) - temperatureold(i, j-1, k, l, m))*(1-exp(-(1/tau12)))); temperaturenew(i, j-1, k, l, m) = temperaturenew(i, j-1, k, l, m) + ((temperatureold(i, j, k, l, m) - temperatureold(i, j-1, k, l, m))*(1-exp(-(1/tau21)))); } if(m == 0) { area = (sqrt(pow((depth-innerdepth)/2, 2) + pow((width-innerwidth)/2, 2))* height) / 0.9 * meterperunit; thermalresistance = sqrt(pow((width/2 - (width-innerwidth)/4)/2, 2)+ pow((depth/2 - (depth-innerdepth)/4)/2, 2))/ 0.9 * meterperunit /(conductivity * area); capacity = specificcapacity *((((depth - innerdepth)/2)*((width + innerwidth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((width - innerwidth)/2)*((depth + innerdepth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; tau21 = capacity2 * thermalresistance; temperaturenew(i, j, k, l, 0) = temperaturenew(i, j, k, l, 0) - ((temperatureold(i, j, k, l, 0) - temperatureold(i, j, k, l, 1))*(1-exp(-(1/tau12)))); temperaturenew(i, j, k, l, 1) = temperaturenew(i, j, k, l, 1) + ((temperatureold(i, j, k, l, 0) - temperatureold(i, j, k, l, 1))*(1-exp(-(1/tau21)))); } else if(m == 1) { area = (sqrt(pow((width-innerwidth)/2, 2) + pow((depth-innerdepth)/2, 2))* height) / 0.9 * meterperunit; thermalresistance = sqrt(pow((depth/2 - (depth-innerdepth)/4)/2, 2)+ pow((width/2 - (width-innerwidth)/4)/2, 2))/ 0.9 * meterperunit /(conductivity * area); capacity = specificcapacity *((((width - innerwidth)/2)*((depth + innerdepth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((depth - innerdepth)/2)*((width + innerwidth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; tau21 = capacity2 * thermalresistance; temperaturenew(i, j, k, l, 1) = temperaturenew(i, j, k, l, 1) - ((temperatureold(i, j, k, l, 1) - temperatureold(i, j, k, l, 2))*(1-exp(-(1/tau12)))); temperaturenew(i, j, k, l, 2) = temperaturenew(i, j, k, l, 2) + ((temperatureold(i, j, k, l, 1) - temperatureold(i, j, k, l, 2))*(1-exp(-(1/tau21)))); } else if(m == 2) { area = (sqrt(pow((depth-innerdepth)/2, 2) + pow((width-innerwidth)/2, 2))* height) / 0.9 * meterperunit; thermalresistance = sqrt(pow((width/2 - (width-innerwidth)/4)/2, 2)+ pow((depth/2 - (depth-innerdepth)/4)/2, 2))/ 0.9 * meterperunit /(conductivity * area); capacity = specificcapacity *((((depth - innerdepth)/2)*((width + innerwidth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((width - innerwidth)/2)*((depth + innerdepth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; tau21 = capacity2 * thermalresistance; temperaturenew(i, j, k, l, 2) = temperaturenew(i, j, k, l, 2) - ((temperatureold(i, j, k, l, 2) - temperatureold(i, j, k, l, 3))*(1-exp(-(1/tau12)))); temperaturenew(i, j, k, l, 3) = temperaturenew(i, j, k, l, 3) + ((temperatureold(i, j, k, l, 2) - temperatureold(i, j, k, l, 3))*(1-exp(-(1/tau21)))); } else if(m == 3) { area = (sqrt(pow((width-innerwidth)/2, 2) + pow((depth-innerdepth)/2, 2))* height) / 0.9 * meterperunit; thermalresistance = sqrt(pow((depth/2 - (depth-innerdepth)/4)/2, 2)+ pow((width/2 - (width-innerwidth)/4)/2, 2))/ 0.9 * meterperunit /(conductivity * area); capacity = specificcapacity *((((width - innerwidth)/2)*((depth + innerdepth)/2)* height)/0.9*meterperunit); capacity2 = specificcapacity *((((depth - innerdepth)/2)*((width + innerwidth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; tau21 = capacity2 * thermalresistance; temperaturenew(i, j, k, l, 3) = temperaturenew(i, j, k, l, 3) - ((temperatureold(i, j, k, l, 3) - temperatureold(i, j, k, l, 0))*(1-exp(-(1/tau12)))); temperaturenew(i, j, k, l, 0) = temperaturenew(i, j, k, l, 0) + ((temperatureold(i, j, k, l, 3) - temperatureold(i, j, k, l, 0))*(1-exp(-(1/tau21)))); } if(j == 0) { if(m % 2 == 0)area = (((depth - innerdepth)/2)*((width + innerwidth)/2)) / 0.9 * meterperunit; if(m % 2 == 1)area = (((width - innerwidth)/2)*((depth + innerdepth)/2)) / 0.9 * meterperunit; thermalresistance = (3 + height/2) / 0.9 * meterperunit /(conductivity * area); if(m % 2 == 0)capacity = specificcapacity *((((depth - innerdepth)/2)* ((width + innerwidth)/2)*height) / 0.9 * meterperunit); if(m % 2 == 1)capacity = specificcapacity *((((width - innerwidth)/2)* ((depth + innerdepth)/2)*height) / 0.9 * meterperunit); tau12 = capacity * thermalresistance; temperaturenew(i, j, k, l, m) = temperaturenew(i, j, k, l, m) - ((temperatureold(i, j, k, l, m) - earthtemp) *(1-exp(-(1/tau12)))); } if(j == sy-1) { if(m % 2 == 0)area = (((depth - innerdepth)/2)*((width + innerwidth)/2)) / 0.9 * meterperunit; if(m % 2 == 1)area = (((width - innerwidth)/2)*((depth + innerdepth)/2)) / 0.9 * meterperunit; thermalresistance = (3 + height/2) / 0.9 * meterperunit /(conductivity * area); if(m % 2 == 0)capacity = specificcapacity *((((depth - innerdepth)/2)* ((width + innerwidth)/2)*height) / 0.9 * meterperunit); if(m % 2 == 1)capacity = specificcapacity *((((width - innerwidth)/2)* ((depth + innerdepth)/2)*height) / 0.9 * meterperunit); tau12 = capacity * thermalresistance; temperaturenew(i, j, k, l, m) = temperaturenew(i, j, k, l, m) - ((temperatureold(i, j, k, l, m) - basementtemp) *(1-exp(-(1/tau12)))); } } } width = cubehole(i, j, k, 0).getWidth(); height = cubehole(i, j, k, 0).getHeight(); depth = cubehole(i, j, k, 0).getDepth(); innerwidth = cubehole(i, j, k, 0).getInnerWidth(); innerdepth = cubehole(i, j, k, 0).getInnerDepth(); if(i > 0) { area = (((depth + innerdepth)/2)*height) / 0.9 * meterperunit; thermalresistance = ((width - innerwidth)/2) / 0.9 * meterperunit / (conductivity * area); capacity = specificcapacity *((((depth + innerdepth)/2)*((width - innerwidth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; temperaturenew(i, j, k, 0, 1) = temperaturenew(i, j, k, 0, 1) - ((temperatureold(i, j, k, 0, 1) - temperatureold(i-1, j, k, 0, 3))*(1-exp(-(1/tau12)))); temperaturenew(i-1, j, k, 0, 3) = temperaturenew(i-1, j, k, 0, 3) + ((temperatureold(i, j, k, 0, 1) - temperatureold(i-1, j, k, 0, 3))*(1-exp(-(1/tau12)))); } if (i == 0 || i == sx-1) { int h=1; if(i == sx-1) h=3; area = (((depth + innerdepth)/2)*height) / 0.9 * meterperunit; thermalresistance = (3 + (width - innerwidth)/4) / 0.9 * meterperunit / (conductivity * area); capacity = specificcapacity *((((depth + innerdepth)/2)*((width - innerwidth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; temperaturenew(i, j, k, 0, h) = temperaturenew(i, j, k, 0, h) - ((temperatureold(i, j, k, 0, h) - earthtemp) *(1-exp(-(1/tau12)))); } if(k > 0) { area = (((width + innerwidth)/2)*height) / 0.9 * meterperunit; thermalresistance = ((depth - innerdepth)/2) / 0.9 * meterperunit / (conductivity * area); capacity = specificcapacity *((((width + innerwidth)/2)*((depth - innerdepth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; temperaturenew(i, j, k, 0, 2) = temperaturenew(i, j, k, 0, 2) - ((temperatureold(i, j, k, 0, 2) - temperatureold(i, j, k-1, 0, 0))*(1-exp(-(1/tau12)))); temperaturenew(i, j, k-1, 0, 0) = temperaturenew(i, j, k-1, 0, 0) + ((temperatureold(i, j, k, 0, 0) - temperatureold(i, j, k-1, 0, 0))*(1-exp(-(1/tau12)))); } if(k == 0 || k == sz-1) { int h=2; if(k == sz-1) h=0; area = (((width + innerwidth)/2)*height) / 0.9 * meterperunit; thermalresistance = (3 + (depth - innerdepth)/4) / 0.9 * meterperunit / (conductivity * area); capacity = specificcapacity *((((width + innerwidth)/2)*((depth - innerdepth)/2)* height)/0.9*meterperunit); tau12 = capacity * thermalresistance; temperaturenew(i, j, k, 0, h) = temperaturenew(i, j, k, 0, h) - ((temperatureold(i, j, k, 0, h) - earthtemp) *(1-exp(-(1/tau12)))); } } } } mergetemperature(); coloring(); // std::cerr << temperaturenew(2, 2, 2, 1, 2) << std::endl; // std::cerr << temperaturenew(2, 2, 2, 2, 2) << std::endl; // std::cerr << temperaturenew(2, 2, 2, 3, 2) << "\n" << std::endl; // std::cerr << " " << temperatureold(2, 2, 2, 1, 2) << std::endl; // std::cerr << temperatureold(2, 2, 2, 2, 1) << " " << temperatureold(2, 2, 2, 2, 2) << " " << temperatureold(2, 2, 2, 2, 3) << std::endl; // std::cerr << " " << temperatureold(2, 2, 2, 3, 2) << "\n" << std::endl; } std::list Temparray::getTriangles(){ std::list triangles; for(int i = 0; i < sx; ++i) { for(int j = 0; j < sy; ++j) { for(int k = 0; k < sz; ++k) { for(int l = 0; l < 6; ++l) { std::list t = cubehole(i, j, k, l).getTriangles(); triangles.splice(triangles.end(), t); } } } } return triangles; } void Temparray::coloring() { calcAverage(); float r[4], g[4], b[4], percent, factor; for(int i=0; i b[m]){ factor = 1 + (1-g[m]) / g[m]; } else { factor = 1 + (1-b[m]) / b[m]; } g[m] *= factor; b[m] *= factor; } else if (percent == 0.5) { r[m] = 0; g[m] = 1; b[m] = 0; } else if (percent > 0.5) { r[m] = (percent-0.5)*2; g[m] = 1-(percent-0.5)*2; b[m] = 0; if(r[m] > g[m]){ factor = 1 + (1-r[m]) / r[m]; } else { factor = 1 + (1-g[m]) / g[m]; } r[m] *= factor; g[m] *= factor; } cubehole(i, j, k, l).setColor(vmml::vec4f(r[0], g[0], b[0], 1.0), vmml::vec4f(r[1], g[1], b[1], 1.0), vmml::vec4f(r[2], g[2], b[2], 1.0), vmml::vec4f(r[3], g[3], b[3], 1.0)); } } } } } } void Temparray::calcweight() { int dist, disth, disth2, posx, posy, posz, posx2, posy2, posz2; maxweight = 0; for(int i=0; i < sx; ++i) { for(int j=0; j < sy; ++j) { for(int k=0; k < sz; ++k) { dist = 0; if (sx%2 == 0) { posx = sx/2.0-1; posx2 = sx/2.0; disth = abs(posx-i); disth2 = abs(posx2-i); if(disth < disth2) dist += disth; else dist += disth2; if (sy%2 == 0) { posy = sy/2.0-1; posy2 = sy/2.0; disth = abs(posy-j); disth2 = abs(posy2-j); if(disth < disth2) dist += disth; else dist += disth2; } else if (sy%2 == 1) { posy = sy/2.0-0.5; dist += abs(posy-j); if (sz%2 == 0) { posz = sz/2.0-1; posz2 = sz/2.0; disth = abs(posz-k); disth2 = abs(posz2-k); if(disth < disth2) dist += disth; else dist += disth2; } else if (sz%2 == 1) { posz = sz/2.0-0.5; dist += abs(posz-k); } } } else if (sx%2 == 1) { posx = sx/2.0-0.5; dist += abs(posx-i); if (sy%2 == 0) { posy = sy/2.0-1; posy2 = sy/2.0; disth = abs(posy-j); disth2 = abs(posy2-j); if(disth < disth2) dist += disth; else dist += disth2; } else if (sy%2 == 1) { posy = sy/2.0-0.5; dist += abs(posy-j); if (sz%2 == 0) { posz = sz/2.0-1; posz2 = sz/2.0; disth = abs(posz-k); disth2 = abs(posz2-k); if(disth < disth2) dist += disth; else dist += disth2; } else if (sz%2 == 1) { posz = sz/2.0-0.5; dist += abs(posz-k); } } } weight(i, j, k) = dist; if (dist > maxweight) maxweight = dist; } } } }