58 #ifdef __TM_MODULE //OOFEM transport module 91 #ifdef __TM_MODULE //OOFEM transport module 113 if ( mode == VM_Total || mode == VM_TotalIntrinsic ) {
124 if ( mode == VM_Total || mode == VM_TotalIntrinsic ) {
171 double conduct = 0.0;
189 if ( !this->
nowarnings.
at(2) && ( conduct < 0.3 || conduct > 5 ) ) {
190 OOFEM_WARNING(
"Weird concrete thermal conductivity %f W/m/K\n", conduct);
202 double capacityConcrete = 0.0;
213 capacityConcrete = 1000 * ms->
GiveCp();
221 if ( !this->
nowarnings.
at(3) && ( capacityConcrete < 500 || capacityConcrete > 2000 ) ) {
222 OOFEM_WARNING(
"Weird concrete heat capacity %f J/kg/K\n", capacityConcrete);
227 return capacityConcrete;
233 double concreteBulkDensity = 0.0;
250 if ( !this->
nowarnings.
at(1) && ( concreteBulkDensity < 1000 || concreteBulkDensity > 4000 ) ) {
251 OOFEM_WARNING(
"Weird concrete density %f kg/m3\n", concreteBulkDensity);
254 concreteBulkDensity *= this->
scaling.
at(1);
256 return concreteBulkDensity;
262 if ( mode == Capacity ) {
264 }
else if ( mode == IntSource ) {
268 double lastEquilibratedTemperature = status->
giveField().
at(1);
270 double krate, EaOverR, val;
276 EaOverR = 1000. * ms->
E_act / 8.314;
278 if ( ms->
icyc > 1 ) {
279 krate = exp( -EaOverR * ( 1. / ( ms->
temp_cur + 273.15 ) - 1. / ( lastEquilibratedTemperature + 273.15 ) ) );
291 val = EaOverR * krate / ( ms->
temp_cur + 273.15 ) / ( ms->
temp_cur + 273.15 );
312 if ( type == IST_HydrationDegree ) {
316 }
else if ( type == IST_Density ) {
320 }
else if ( type == IST_ThermalConductivityIsotropic ) {
324 }
else if ( type == IST_HeatCapacity ) {
328 }
else if ( type == IST_AverageTemperature ) {
332 }
else if ( type == IST_YoungModulusVirginPaste ) {
336 }
else if ( type == IST_PoissonRatioVirginPaste ) {
340 }
else if ( type == IST_YoungModulusConcrete ) {
344 }
else if ( type == IST_PoissonRatioConcrete ) {
368 gp->setMaterialStatus( ms, this->
giveNumber() );
393 element->
giveIPValue(vecTemperature, gp, IST_Temperature, tStep);
413 if ( result !=
IRRT_OK )
return result;
452 OOFEM_ERROR(
"Use function CemhydMat :: initMaterial instead");
538 printf(
"Constructor of CemhydMatStatus on GP %p, withMicrostructure %d, copy from CemhydMatStatus %p\n", gp, withMicrostructure, CemStat);
542 if ( withMicrostructure ) {
548 for ( k = 0; k <
SYSIZE; k++ ) {
549 for ( j = 0; j <
SYSIZE; j++ ) {
550 for ( i = 0; i <
SYSIZE; i++ ) {
553 mic [ i ] [ j ] [ k ] =
micorig [ i ] [ j ] [ k ];
827 phase =
new long int [ 51 ];
849 for (
int i = 0; i <=
EMPTYP; i++ ) {
854 for (
int i = 0; i <=
HDCSH; i++ ) {
956 printf(
"Deallocating headant\n");
966 printf(
"Deallocating arrays\n");
984 for (
int x = 0; x <
SYSIZE; x++ ) {
986 for (
int y = 0; y <
SYSIZE; y++ ) {
988 if (
mic [ x ] [ y ] == NULL ) {
989 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
998 for (
int x = 0; x <
SYSIZE; x++ ) {
999 for (
int y = 0; y <
SYSIZE; y++ ) {
1000 delete []
mic [ x ] [ y ];
1003 delete []
mic [ x ];
1013 for (
int x = 0; x <
SYSIZE; x++ ) {
1015 for (
int y = 0; y <
SYSIZE; y++ ) {
1017 if (
mic [ x ] [ y ] == NULL ) {
1018 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1027 if (
mic != NULL ) {
1028 for (
int x = 0; x <
SYSIZE; x++ ) {
1029 for (
int y = 0; y <
SYSIZE; y++ ) {
1030 delete []
mic [ x ] [ y ];
1033 delete []
mic [ x ];
1043 for (
int x = 0; x <
SYSIZE; x++ ) {
1045 for (
int y = 0; y <
SYSIZE; y++ ) {
1047 if (
mic [ x ] [ y ] == NULL ) {
1048 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1057 if (
mic != NULL ) {
1058 for (
int x = 0; x <
SYSIZE; x++ ) {
1059 for (
int y = 0; y <
SYSIZE; y++ ) {
1060 delete []
mic [ x ] [ y ];
1063 delete []
mic [ x ];
1073 for (
int x = 0; x <
SYSIZE; x++ ) {
1075 for (
int y = 0; y <
SYSIZE; y++ ) {
1076 mic [ x ] [ y ] =
new short int [
SYSIZE ];
1077 if (
mic [ x ] [ y ] == NULL ) {
1078 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1087 if (
mic != NULL ) {
1088 for (
int x = 0; x <
SYSIZE; x++ ) {
1089 for (
int y = 0; y <
SYSIZE; y++ ) {
1090 delete []
mic [ x ] [ y ];
1093 delete []
mic [ x ];
1103 for (
int x = 0; x <
SYSIZE; x++ ) {
1105 for (
int y = 0; y <
SYSIZE; y++ ) {
1106 mic [ x ] [ y ] =
new double [
SYSIZE ];
1107 if (
mic [ x ] [ y ] == NULL ) {
1108 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1117 if (
mic != NULL ) {
1118 for (
int x = 0; x <
SYSIZE; x++ ) {
1119 for (
int y = 0; y <
SYSIZE; y++ ) {
1120 delete []
mic [ x ] [ y ];
1123 delete []
mic [ x ];
1136 XMLHandle docHandle = XMLHandle(xmlFile);
1137 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1138 if ( elemSelected == NULL ) {
1139 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1143 sprintf(key,
"key%d", position);
1144 success = elemSelected->QueryIntAttribute(key, & val);
1145 if ( success != XML_SUCCESS ) {
1146 printf(
"Cannot read int value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1155 val =
static_cast< long int >(temp);
1161 XMLHandle docHandle = XMLHandle(xmlFile);
1162 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1163 if ( elemSelected == NULL ) {
1164 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1168 success = elemSelected->QueryIntAttribute(key, & val);
1169 if ( success != XML_SUCCESS ) {
1170 printf(
"Cannot read int value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1180 XMLHandle docHandle = XMLHandle(xmlFile);
1181 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1182 if ( elemSelected == NULL ) {
1183 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1187 sprintf(key,
"key%d", position);
1188 success = elemSelected->QueryDoubleAttribute(key, & val);
1189 if ( success != XML_SUCCESS ) {
1190 printf(
"Cannot read double value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1198 XMLHandle docHandle = XMLHandle(xmlFile);
1199 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1200 if ( elemSelected == NULL ) {
1201 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1205 success = elemSelected->QueryDoubleAttribute(key, & val);
1206 if ( success != XML_SUCCESS ) {
1207 printf(
"Cannot read double value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1216 XMLHandle docHandle = XMLHandle(xmlFile);
1218 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1219 if ( elemSelected == NULL ) {
1220 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1224 sprintf(key,
"key%d", position);
1228 const char *cstr = elemSelected->Attribute(key);
1230 success = XML_SUCCESS;
1232 success = XML_NO_ATTRIBUTE;
1234 if ( success != XML_SUCCESS ) {
1235 printf(
"Cannot read string value or key %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1238 str1 = std :: string(cstr);
1240 strcpy( chars, str1.c_str() );
1252 F =
new cmlfile(inp);
1255 F->set_label_string(0,
"Rand_seed_num");
1256 F->set_label_string(1,
"Input_img_file");
1257 F->set_label_string(2,
"Input_id_file");
1258 F->set_label_string(3,
"Saturated_sealed");
1259 F->set_label_string(4,
"Induction_time");
1260 F->set_label_string(5,
"Ea_cement");
1261 F->set_label_string(6,
"Ea_pozz");
1262 F->set_label_string(7,
"Ea_slag");
1263 F->set_label_string(8,
"Beta");
1264 F->set_label_string(9,
"Mass_SCM_FA_CA_inert_frac");
1265 F->set_label_string(10,
"Mass_cem");
1266 F->set_label_string(11,
"Cp_SCM_FA_CA_inert");
1267 F->set_label_string(12,
"Cp_cem");
1268 F->set_label_string(13,
"Given_microstructure");
1269 F->set_label_string(14,
"Output_initial_microstructure");
1270 F->set_label_string(15,
"Output_initial_microstructure_img_file");
1271 F->set_label_string(16,
"Output_initial_microstructure_id_file");
1272 F->set_label_string(17,
"Cycle_freq_perc_pore");
1273 F->set_label_string(18,
"Cycle_freq_perc_sol");
1274 F->set_label_string(19,
"Total_sodium");
1275 F->set_label_string(20,
"Total_potassium");
1276 F->set_label_string(21,
"Readily_soluble_sodium");
1277 F->set_label_string(22,
"Readily_soluble_potassium");
1278 F->set_label_string(23,
"Diffusion_steps_per_cycle");
1279 F->set_label_string(24,
"CH_nucleation_probability");
1280 F->set_label_string(25,
"CH_scale_factor");
1281 F->set_label_string(26,
"Gypsum_nucleation_probability");
1282 F->set_label_string(27,
"Gypsum_scale_factor");
1283 F->set_label_string(28,
"C3AH6_nucleation_probability");
1284 F->set_label_string(29,
"C3AH6_scale_factor");
1285 F->set_label_string(30,
"FH3_nucleation_probability");
1286 F->set_label_string(31,
"FH3_scale_factor");
1287 F->set_label_string(32,
"Microstructure_size");
1288 F->set_label_string(33,
"Adiabatic_conditions");
1289 F->set_label_string(34,
"Vol_cement_clinker_gypsum");
1290 F->set_label_string(35,
"Vol_cement_SCM");
1291 F->set_label_string(36,
"Vol_water");
1292 F->set_label_string(37,
"Vol_FA");
1293 F->set_label_string(38,
"Vol_CA");
1294 F->set_label_string(39,
"Vol_inert_filler");
1295 F->set_label_string(40,
"Vol_entrained_entrapped_air");
1296 F->set_label_string(41,
"Grain_average_FA");
1297 F->set_label_string(42,
"Grain_average_CA");
1298 F->set_label_string(43,
"ITZ_thickness");
1299 F->set_label_string(44,
"ITZ_Young_red");
1300 F->set_label_string(45,
"Young_SCM");
1301 F->set_label_string(46,
"Poisson_SCM");
1302 F->set_label_string(47,
"Young_FA");
1303 F->set_label_string(48,
"Poisson_FA");
1304 F->set_label_string(49,
"Young_CA");
1305 F->set_label_string(50,
"Poisson_CA");
1306 F->set_label_string(51,
"Young_inert");
1307 F->set_label_string(52,
"Poisson_inert");
1308 F->set_label_string(53,
"Calculate_elastic_homogenization");
1363 F->set_section_string(0,
"CEMHYD_generate_particles");
1365 F->check_requirements();
1367 if ( F->error_in_requirements() ) {
1368 printf(
"Cemhyd input file %s is not complete (file %s, line %d)\n", inp, __FILE__, __LINE__);
1372 F->get_value(0, (
long & )
iseed);
1377 int errorId =
xmlFile->LoadFile(inp);
1378 if ( errorId != XML_SUCCESS ) {
1379 printf(
"\nError reading XML file %s or nonletter symbols used, error id = %d\n", inp, errorId);
1394 F->get_value(32, (
long & )
SYSSIZE);
1396 if ( SYSSIZE < 10 ) {
1397 printf(
"Can not run small microstructure %d (< 10 voxels a side), file %s, line %d)\n", SYSSIZE, __FILE__, __LINE__);
1428 F->get_value(13, (
long & )read_micr);
1431 if ( !read_micr && generateMicrostructure == 1 ) {
1437 printf(
"MONOPHASE microstructure created\n");
1461 if ( ( * idum <= 0 ) || (
iy == 0 ) ) {
1462 * idum = ( -* idum > * idum ) ? -* idum : * idum;
1463 for ( j =
NTAB + 7; j >= 0; j-- ) {
1465 * idum =
IA * ( * idum - k *
IQ ) -
IR * k;
1479 * idum =
IA * ( * idum - k *
IQ ) -
IR * k;
1484 j = ( int ) (
iy *
NDIV );
1501 printf(
"Enter thickness of aggregate to place (an even integer) \n");
1503 F->get_next_line_in_section(0, (
long & )
aggsize);
1511 printf(
"%d\n", aggsize);
1520 for ( ix = agglo; ix <= agghi; ix++ ) {
1521 for ( iy = 1; iy <=
SYSSIZE; iy++ ) {
1522 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
1543 int nofits, xp, yp, zp, i, j, k;
1544 float dist, xdist, ydist, zdist, ftmp;
1548 for ( i = xin - radd; ( ( i <= xin + radd ) && ( nofits == 0 ) ); i++ ) {
1557 ftmp = ( float ) ( i - xin );
1558 xdist = ftmp * ftmp;
1559 for ( j = yin - radd; ( ( j <= yin + radd ) && ( nofits == 0 ) ); j++ ) {
1568 ftmp = ( float ) ( j - yin );
1569 ydist = ftmp * ftmp;
1570 for ( k = zin - radd; ( ( k <= zin + radd ) && ( nofits == 0 ) ); k++ ) {
1579 ftmp = ( float ) ( k - zin );
1580 zdist = ftmp * ftmp;
1583 dist = sqrt(xdist + ydist + zdist);
1584 if ( ( dist - 0.5 ) <= (
float ) radd ) {
1587 cement [ xp ] [ yp ] [ zp ] = phasein;
1588 cemreal [ xp ] [ yp ] [ zp ] = phase2;
1591 else if ( ( wflg == 1 ) && (
cement [ xp ] [ yp ] [ zp ] !=
POROSITY ) ) {
1597 if ( ( wflg == 1 ) && ( ( fabs( xp - ( (
float ) (
SYSSIZE + 1 ) / 2.0 ) ) ) < ( (
float )
aggsize / 2.0 ) ) ) {
1620 int count, x, y, z, radius, ig, tries, phnow;
1622 float testgyp, typegyp;
1625 for ( ig = 0; ig < numgen; ig++ ) {
1626 phnow = pheach [ ig ];
1627 radius = sizeeach [ ig ];
1629 for ( jg = 1; jg <= numeach [ ig ]; jg++ ) {
1648 printf(
"Could not place sphere %d after %ld random attempts \n",
npart,
MAXTRIES);
1649 printf(
"Skipping this microstructure parameters\n");
1650 for ( i = 1; i <=
npart; i++ ) {
1656 }
while ( count != 0 );
1661 printf(
"Too many spheres being generated \n");
1662 printf(
"User needs to increase value of NPARTC at top of C-code\n");
1663 printf(
"Skipping this microstructure parameters\n");
1718 for ( i = 1; i <=
npart; i++ ) {
1733 int *sphrad, *sphase;
1740 sphnum =
new long int [
NUMSIZES ];
1745 printf(
"Enter number of different size spheres to use(max. is %d) \n",
NUMSIZES);
1748 F->get_next_line_in_section(0, (
long & )numsize);
1755 printf(
"%d \n", numsize);
1757 }
while ( ( numsize >
NUMSIZES ) || ( numsize < 0 ) );
1761 printf(
"Enter dispersion factor (separation distance in pixels) for spheres (0-2) \n");
1762 printf(
"0 corresponds to totally random placement \n");
1765 F->get_next_line_in_section(0, (
long & )
dispdist);
1772 printf(
"%d \n", dispdist);
1778 printf(
"Enter probability for gypsum particles on a random particle basis (0.0-1.0) \n");
1781 F->get_next_line_in_section(0,
probgyp);
1794 printf(
"Enter probabilities for hemihydrate and anhydrite forms of gypsum (0.0-1.0) \n");
1797 F->get_next_line_in_section(0,
probhem);
1798 F->get_next_line_in_section(0,
probanh);
1810 if ( ( numsize > 0 ) && ( numsize < (
NUMSIZES + 1 ) ) ) {
1812 printf(
"Enter number, radius, and phase ID for each sphere class (largest radius 1st) \n");
1813 printf(
"Phases are %d- Cement and (random) calcium sulfate, %d- C2S, %d- Gypsum, %d- hemihydrate %d- anhydrite %d- Pozzolanic, %d- Inert, %d- Slag, %d- CaCO3 %d- Fly Ash \n",
CEMID,
C2SID,
GYPID,
HEMIHYDRATE,
ANHYDRITE,
POZZID,
INERTID,
SLAGID,
CACO3,
FLYASH);
1816 for ( isph = 0; isph < numsize; isph++ ) {
1818 printf(
"Enter number of spheres of class %d \n", isph + 1);
1821 F->get_next_line_in_section(0, inval1);
1829 printf(
"%ld \n", inval1);
1831 sphnum [ isph ] = inval1;
1835 printf(
"Enter radius of spheres of class %d \n", isph + 1);
1836 printf(
"(Integer <=%d please) \n",
SYSSIZE / 3);
1839 F->get_next_line_in_section(0, (
long & )inval);
1845 if ( inval > (
SYSSIZE / 3 ) ) {
1846 printf(
"Given radius %d exceeded maximum radius of %d, terminating\n", inval,
SYSSIZE / 3);
1851 printf(
"%d \n", inval);
1855 sphrad [ isph ] = inval;
1858 printf(
"Enter phase of spheres of class %d \n", isph + 1);
1861 F->get_next_line_in_section(0, (
long & )inval);
1868 printf(
"%d \n", inval);
1872 sphase [ isph ] = inval;
1873 if ( inval ==
CEMID ) {
1882 if (
gsphere(numsize, sphnum, sphrad, sphase) == 1 ) {
1905 int xp, yp, zp, i, j, k;
1906 float dist, xdist, ydist, zdist, ftmp;
1909 for ( i = xin - radd; ( i <= xin + radd ); i++ ) {
1918 ftmp = ( float ) ( i - xin );
1919 xdist = ftmp * ftmp;
1920 for ( j = yin - radd; ( j <= yin + radd ); j++ ) {
1929 ftmp = ( float ) ( j - yin );
1930 ydist = ftmp * ftmp;
1931 for ( k = zin - radd; ( k <= zin + radd ); k++ ) {
1940 ftmp = ( float ) ( k - zin );
1941 zdist = ftmp * ftmp;
1944 dist = sqrt(xdist + ydist + zdist);
1945 if ( ( dist - 0.5 ) <= (
float ) radd ) {
1947 cement [ xp ] [ yp ] [ zp ] = phasein;
1948 cemreal [ xp ] [ yp ] [ zp ] = phase2;
1963 int nofits, xp, yp, zp, i, j, k;
1964 float dist, xdist, ydist, zdist, ftmp;
1969 for ( i = xin - radd; ( ( i <= xin + radd ) && ( nofits == 0 ) ); i++ ) {
1978 ftmp = ( float ) ( i - xin );
1979 xdist = ftmp * ftmp;
1980 for ( j = yin - radd; ( ( j <= yin + radd ) && ( nofits == 0 ) ); j++ ) {
1989 ftmp = ( float ) ( j - yin );
1990 ydist = ftmp * ftmp;
1991 for ( k = zin - radd; ( ( k <= zin + radd ) && ( nofits == 0 ) ); k++ ) {
2000 ftmp = ( float ) ( k - zin );
2001 zdist = ftmp * ftmp;
2004 dist = sqrt(xdist + ydist + zdist);
2005 if ( ( dist - 0.5 ) <= (
float ) radd ) {
2008 nofits =
cement [ xp ] [ yp ] [ zp ];
2013 if ( ( fabs( xp - ( (
float ) (
SYSSIZE + 1 ) / 2.0 ) ) ) < ( (
float )
aggsize / 2.0 ) ) {
2030 int partdo, numfloc;
2032 int nleft = 0, ckall;
2033 int xm, ym, zm, moveran;
2034 int xp, yp, zp, rp, clushit, valkeep;
2037 struct cluster *parttmp, *partpoint, *partkeep = NULL;
2039 cluspart =
new int [
NPARTC ];
2043 for ( iclus = 1; iclus <=
npart; iclus++ ) {
2044 cluspart [ iclus ] = iclus;
2049 printf(
"Enter number of flocs desired at end of routine (>0) \n");
2052 F->get_next_line_in_section(0, (
long & )numfloc);
2059 printf(
"%d\n", numfloc);
2061 }
while ( numfloc <= 0 );
2063 while ( nstart > numfloc ) {
2067 for ( iclus = 1; iclus <=
npart; iclus++ ) {
2068 if (
clust [ iclus ] == NULL ) {
2073 moveran = ( int ) ( 6. *
ran1(
seed) );
2074 switch ( moveran ) {
2098 partpoint =
clust [ iclus ];
2099 while ( partpoint != NULL ) {
2110 partpoint =
clust [ iclus ];
2111 while ( ( partpoint != NULL ) && ( ckall == 0 ) ) {
2112 xp = partpoint->
x + xm;
2113 yp = partpoint->
y + ym;
2114 zp = partpoint->
z + zm;
2116 ckall =
chkfloc(xp, yp, zp, rp);
2122 partpoint =
clust [ iclus ];
2123 while ( partpoint != NULL ) {
2124 xp = partpoint->
x + xm;
2125 yp = partpoint->
y + ym;
2126 zp = partpoint->
z + zm;
2129 partdo = partpoint->
partid;
2130 drawfloc(xp, yp, zp, rp, partdo +
CEM - 1, valkeep);
2140 partpoint =
clust [ iclus ];
2142 while ( partpoint != NULL ) {
2148 partdo = partpoint->
partid;
2149 drawfloc(xp, yp, zp, rp, partdo +
CEM - 1, valkeep);
2150 partkeep = partpoint;
2155 if ( ckall !=
AGG ) {
2156 clushit = cluspart [ ckall -
CEM + 1 ];
2158 parttmp =
clust [ clushit ];
2161 while ( parttmp != NULL ) {
2162 cluspart [ parttmp->
partid ] = iclus;
2169 clust [ clushit ] = NULL;
2176 printf(
"Number left was %d but number of clusters is %d \n", nleft, nstart);
2191 long int npor, nc2s, ngyp, ncem, nagg, npozz, ninert, nflyash, nanh, nhem, ncaco3, nslag;
2209 for ( i = 1; i <=
SYSSIZE; i++ ) {
2210 for ( j = 1; j <=
SYSSIZE; j++ ) {
2211 for ( k = 1; k <=
SYSSIZE; k++ ) {
2212 valph =
cemreal [ i ] [ j ] [ k ];
2215 }
else if ( valph ==
CEMID ) {
2217 }
else if ( valph ==
C2SID ) {
2219 }
else if ( valph ==
GYPID ) {
2225 }
else if ( valph ==
AGG ) {
2227 }
else if ( valph ==
POZZID ) {
2229 }
else if ( valph ==
SLAGID ) {
2231 }
else if ( valph ==
INERTID ) {
2233 }
else if ( valph ==
FLYASH ) {
2235 }
else if ( valph ==
CACO3 ) {
2243 printf(
"\n Phase counts are: \n");
2244 printf(
"Porosity= %ld \n", npor);
2245 printf(
"Cement= %ld \n", ncem);
2246 printf(
"C2S= %ld \n", nc2s);
2247 printf(
"Gypsum= %ld \n", ngyp);
2248 printf(
"Anhydrite= %ld \n", nanh);
2249 printf(
"Hemihydrate= %ld \n", nhem);
2250 printf(
"Pozzolan= %ld \n", npozz);
2251 printf(
"Inert= %ld \n", ninert);
2252 printf(
"Slag= %ld \n", nslag);
2253 printf(
"CaCO3= %ld \n", ncaco3);
2254 printf(
"Fly Ash= %ld \n", nflyash);
2255 printf(
"Aggregate= %ld \n", nagg);
2265 int phase [ 40 ], ptot;
2266 int icnt, ixlo, ixhi,
iy, iz, phid, idist;
2270 aggfile = fopen(
"agglist.out",
"w");
2271 printf(
"Distance Porosity Cement C2S Gypsum Anhydrite Hemihydrate Pozzolan Inert Slag CaCO3 Fly Ash\n");
2272 fprintf(aggfile,
"Distance Porosity Cement C2S Gypsum Anhydrite Hemihydrate Pozzolan Inert Slag CaCO3 Fly Ash\n");
2282 for ( icnt = 0; icnt < 39; icnt++ ) {
2289 for ( iy = 1; iy <=
SYSSIZE; iy++ ) {
2290 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2294 phase [ phid ] += 1;
2300 phase [ phid ] += 1;
2306 printf(
"%d %d %d %d %d %d %d %d %d %d %d %d\n", idist, phase [ 0 ], phase [
CEMID ], phase [
C2SID ], phase [
GYPID ], phase [
ANHYDRITE ], phase [
HEMIHYDRATE ], phase [
POZZID ], phase [
INERTID ], phase [
SLAGID ], phase [
CACO3 ], phase [
FLYASH ]);
2307 fprintf(aggfile,
"%d %d %d %d %d %d %d %d %d %d %d %d\n", idist, phase [ 0 ], phase [ CEMID ], phase [ C2SID ], phase [ GYPID ], phase [ ANHYDRITE ], phase [ HEMIHYDRATE ], phase [ POZZID ], phase [ INERTID ], phase [ SLAGID ], phase [ CACO3 ], phase [ FLYASH ]);
2321 long int inew, ntop, nthrough, ncur, nnew, ntot;
2322 int i, j, k, nmatx [ 29000 ], nmaty [ 29000 ], nmatz [ 29000 ];
2323 int xcn, ycn, zcn, npix, x1, y1, z1, igood, nnewx [ 29000 ], nnewy [ 29000 ], nnewz [ 29000 ];
2327 printf(
"Enter phase to analyze 0) pores 1) Cement \n");
2329 F->get_next_line_in_section(0, (
long & )npix);
2335 printf(
"%d \n", npix);
2336 }
while ( ( npix != 0 ) && ( npix != 1 ) );
2346 for ( i = 1; i <=
SYSSIZE; i++ ) {
2347 for ( j = 1; j <=
SYSSIZE; j++ ) {
2351 if ( ( (
cement [ i ] [ j ] [ k ] == npix ) && ( (
cement [ i ] [ j ] [
SYSSIZE ] == npix ) ||
2367 for ( inew = 1; inew <= ncur; inew++ ) {
2368 xcn = nmatx [ inew ];
2369 ycn = nmaty [ inew ];
2370 zcn = nmatz [ inew ];
2373 for ( jnew = 1; jnew <= 6; jnew++ ) {
2382 }
else if ( jnew == 2 ) {
2387 }
else if ( jnew == 3 ) {
2392 }
else if ( jnew == 4 ) {
2397 }
else if ( jnew == 5 ) {
2399 }
else if ( jnew == 6 ) {
2404 if ( ( z1 >= 1 ) && ( z1 <=
SYSSIZE ) ) {
2405 if ( (
cement [ x1 ] [ y1 ] [ z1 ] == npix ) || ( (
cement [ x1 ] [ y1 ] [ z1 ] >=
CEM ) &&
2406 (
cement [ x1 ] [ y1 ] [ z1 ] <
BURNTG ) && ( npix == 1 ) ) ) {
2410 if ( nnew >= 29000 ) {
2411 printf(
"error in size of nnew \n");
2414 nnewx [ nnew ] = x1;
2415 nnewy [ nnew ] = y1;
2416 nnewz [ nnew ] = z1;
2429 for ( icur = 1; icur <= ncur; icur++ ) {
2430 nmatx [ icur ] = nnewx [ icur ];
2431 nmaty [ icur ] = nnewy [ icur ];
2432 nmatz [ icur ] = nnewz [ icur ];
2435 }
while ( nnew > 0 );
2445 printf(
"Phase ID= %d \n", npix);
2446 printf(
"Number accessible from top= %ld \n", ntop);
2447 printf(
"Number contained in through pathways= %ld \n", nthrough);
2450 for ( i = 1; i <=
SYSSIZE; i++ ) {
2451 for ( j = 1; j <=
SYSSIZE; j++ ) {
2452 for ( k = 1; k <=
SYSSIZE; k++ ) {
2467 FILE *outfile, *partfile;
2468 char filen [ 80 ], filepart [ 80 ];
2469 int ix,
iy, iz, valout;
2472 printf(
"Enter name of file to save microstructure to \n");
2475 F->get_next_line_in_section(0, filen);
2481 printf(
"%s\n", filen);
2483 outfile = fopen(filen,
"w");
2486 printf(
"Enter name of file to save particle IDs to \n");
2489 F->get_next_line_in_section(0, filepart);
2496 printf(
"%s\n", filepart);
2499 partfile = fopen(filepart,
"w");
2501 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2502 for ( iy = 1; iy <=
SYSSIZE; iy++ ) {
2503 for ( ix = 1; ix <=
SYSSIZE; ix++ ) {
2505 fprintf(outfile,
"%1d\n", valout);
2506 valout =
cement [ ix ] [
iy ] [ iz ];
2511 fprintf(partfile,
"%d\n", valout);
2591 for ( ig = 1; ig <=
SYSSIZE; ig++ ) {
2592 for ( jg = 1; jg <=
SYSSIZE; jg++ ) {
2593 for ( kg = 1; kg <=
SYSSIZE; kg++ ) {
2603 printf(
" \n Input User Choice \n");
2604 printf(
"1) Exit \n");
2605 printf(
"2) Add spherical particles (cement,gypsum, pozzolans, etc.) to microstructure \n");
2606 printf(
"3) Flocculate system by reducing number of particle clusters \n");
2607 printf(
"4) Measure global phase fractions \n");
2608 printf(
"5) Add an aggregate to the microstructure \n");
2609 printf(
"6) Measure single phase connectivity (pores or solids) \n");
2610 printf(
"7) Measure phase fractions vs. distance from aggregate surface \n");
2611 printf(
"8) Output current microstructure to file \n");
2614 F->get_next_line_in_section(0, (
long & )userc);
2649 printf(
"No aggregate present. \n");
2659 }
while ( userc != 1 );
2662 for ( ig = 0; ig <
SYSSIZE; ig++ ) {
2663 for ( jg = 0; jg <
SYSSIZE; jg++ ) {
2664 for ( kg = 0; kg <
SYSSIZE; kg++ ) {
2665 micpart [ ig ] [ jg ] [ kg ] =
cement [ ig + 1 ] [ jg + 1 ] [ kg + 1 ];
2686 int icirc, xval, yval, zval;
2692 for ( xval = ( -size ); xval <= size; xval++ ) {
2693 xtmp = ( float ) ( xval * xval );
2694 for ( yval = ( -size ); yval <= size; yval++ ) {
2695 ytmp = ( float ) ( yval * yval );
2696 for ( zval = ( -size ); zval <= size; zval++ ) {
2697 dist = sqrt( xtmp + ytmp + (
float ) ( zval * zval ) );
2698 if ( dist <= ( (
float ) size + 0.5 ) ) {
2701 printf(
"Too many elements in sphere \n");
2702 printf(
"Must change value of MAXSPH parameter \n");
2703 printf(
"Currently set at %ld \n",
MAXSPH);
2707 xsph [ icirc ] = xval;
2708 ysph [ icirc ] = yval;
2709 zsph [ icirc ] = zval;
2724 long int npore,
nsolid [ 37 ];
2728 for ( ix = 1; ix < 37; ix++ ) {
2733 for ( ix = 1; ix <=
SYSSIZE; ix++ ) {
2734 for ( iy = 1; iy <=
SYSSIZE; iy++ ) {
2735 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2736 if (
mask [ ix ] [ iy ] [ iz ] == 0 ) {
2739 nsolid [
mask [ ix ] [
iy ] [ iz ] ] += 1;
2745 printf(
"Pores are: %ld \n", npore);
2746 printf(
"Solids are: %ld %ld %ld %ld %ld %ld\n", nsolid [ 1 ], nsolid [ 2 ],
2747 nsolid [ 3 ], nsolid [ 4 ], nsolid [ 5 ], nsolid [ 6 ]);
2756 int npix, ix1, iy1, iz1;
2767 if (
mask [ ix1 ] [ yin ] [ zin ] == 0 ) {
2776 if (
mask [ ix1 ] [ yin ] [ zin ] == 0 ) {
2785 if (
mask [ xin ] [ iy1 ] [ zin ] == 0 ) {
2794 if (
mask [ xin ] [ iy1 ] [ zin ] == 0 ) {
2803 if (
mask [ xin ] [ yin ] [ iz1 ] == 0 ) {
2812 if (
mask [ xin ] [ yin ] [ iz1 ] == 0 ) {
2825 long int porc, surfc;
2831 for ( ix = 1; ix <=
SYSSIZE; ix++ ) {
2832 for ( iy = 1; iy <=
SYSSIZE; iy++ ) {
2833 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2834 if (
mask [ ix ] [ iy ] [ iz ] == phin ) {
2842 printf(
"Phase area count is %ld \n", porc);
2843 printf(
"Phase surface count is %ld \n", surfc);
2844 rhval = ( float ) porc * 6. / ( 4. * (
float ) surfc );
2845 printf(
"Hydraulic radius is %f \n", rhval);
2859 for ( ic = 1; ic <=
nsph; ic++ ) {
2860 xc = xp +
xsph [ ic ];
2861 yc = yp +
ysph [ ic ];
2862 zc = zp +
zsph [ ic ];
2882 if ( ( xc != xp ) || ( yc != yp ) || ( zc != zp ) ) {
2883 if ( (
mask [ xc ] [ yc ] [ zc ] == phin ) || (
mask [ xc ] [ yc ] [ zc ] == 0 ) ) {
2898 int count, xl, yl, zl;
2902 for ( xl = 1; xl <=
SYSSIZE; xl++ ) {
2903 for ( yl = 1; yl <=
SYSSIZE; yl++ ) {
2904 for ( zl = 1; zl <=
SYSSIZE; zl++ ) {
2908 if (
mask [ xl ] [ yl ] [ zl ] == ph1 ) {
2909 count =
countem(xl, yl, zl, 0);
2914 if (
mask [ xl ] [ yl ] [ zl ] == ph2 ) {
2915 count =
countem(xl, yl, zl, ph2);
2918 if ( ( count < 0 ) || ( count >=
nsph ) ) {
2919 printf(
"Error count is %d \n", count);
2920 printf(
"xl %d yl %d zl %d \n", xl, yl, zl);
2925 if ( ( count >= 0 ) && (
mask [ xl ] [ yl ] [ zl ] == ph1 ) ) {
2932 if ( ( count >= 0 ) && (
mask [ xl ] [ yl ] [ zl ] == ph2 ) ) {
2950 int xd, yd, zd, curvval;
2953 for ( xd = 1; xd <=
SYSSIZE; xd++ ) {
2954 for ( yd = 1; yd <=
SYSSIZE; yd++ ) {
2955 for ( zd = 1; zd <=
SYSSIZE; zd++ ) {
2956 curvval =
curvature [ xd ] [ yd ] [ zd ];
2958 if (
mask [ xd ] [ yd ] [ zd ] == ph2 ) {
2959 nair [ curvval ] += 1;
2960 }
else if (
mask [ xd ] [ yd ] [ zd ] == ph1 ) {
2975 int valfound, i, stop;
2980 valfound =
nsph - 1;
2983 for ( i = (
nsph - 1 ); ( ( i >= 0 ) && ( stop == 0 ) ); i-- ) {
2985 if ( nsofar > nsearch ) {
2991 return ( valfound );
3002 int valfound, i, stop;
3010 for ( i = 0; ( ( i <
nsph ) && ( stop == 0 ) ); i++ ) {
3011 nsofar +=
nair [ i ];
3012 if ( nsofar > nsearch ) {
3018 return ( valfound );
3027 int count1, count2, ntot, countc, i, xp, yp, zp;
3028 int cmin, cmax, cfg;
3030 long int nsolc, nairc, nsum, nsolm, nairm, nst1, nst2, next1, next2;
3031 float pck, plsol, plair;
3039 for ( i =
nsph; i > count1; i-- ) {
3040 if ( (
nsolid [ i ] > 0 ) && ( cfg == 0 ) ) {
3049 plsol = ( float ) ( ntomove - nsum ) / ( float )
nsolid [ count1 ];
3050 next1 = ntomove - nsum;
3051 nst1 =
nsolid [ count1 ];
3057 for ( i = 0; i < count2; i++ ) {
3058 if ( (
nair [ i ] > 0 ) && ( cfg == 0 ) ) {
3067 plair = ( float ) ( ntomove - nsum ) / ( float )
nair [ count2 ];
3068 next2 = ntomove - nsum;
3069 nst2 =
nair [ count2 ];
3073 if ( cmin >= cmax ) {
3075 printf(
"Stopping - at equilibrium \n");
3076 printf(
"cmin- %d cmax- %d \n", cmin, cmax);
3088 for ( xp = 1; xp <=
SYSSIZE; xp++ ) {
3089 for ( yp = 1; yp <=
SYSSIZE; yp++ ) {
3090 for ( zp = 1; zp <=
SYSSIZE; zp++ ) {
3091 countc =
curvature [ xp ] [ yp ] [ zp ];
3093 if (
mask [ xp ] [ yp ] [ zp ] == ph1 ) {
3094 if ( countc > count1 ) {
3096 mask [ xp ] [ yp ] [ zp ] = ph2;
3100 nair [ countc ] += 1;
3105 if ( countc == count1 ) {
3109 if ( ( pck < 0 ) || ( pck > 1.0 ) ) {
3113 if ( ( ( pck < plsol ) && ( nsolc < next1 ) ) || ( ( nst1 - nsolm ) < ( next1 - nsolc ) ) ) {
3116 mask [ xp ] [ yp ] [ zp ] = ph2;
3120 nair [ count1 ] += 1;
3127 else if (
mask [ xp ] [ yp ] [ zp ] == ph2 ) {
3128 if ( countc < count2 ) {
3130 mask [ xp ] [ yp ] [ zp ] = ph1;
3133 nair [ countc ] -= 1;
3137 if ( countc == count2 ) {
3140 if ( ( pck < 0 ) || ( pck > 1.0 ) ) {
3144 if ( ( ( pck < plair ) && ( nairc < next2 ) ) || ( ( nst2 - nairm ) < ( next2 - nairc ) ) ) {
3147 mask [ xp ] [ yp ] [ zp ] = ph1;
3150 nair [ count2 ] -= 1;
3164 printf(
"ntot is %d \n", ntot);
3173 int natonce, i, rade, j, rflag;
3175 long int curvsum1, curvsum2, pixsum1, pixsum2;
3176 float rhnow, avecurv1, avecurv2;
3179 for ( i = 0; i <= 499; i++ ) {
3190 printf(
"nsph is %d \n",
nsph);
3199 while ( ( rhnow < rhtarget ) && ( i <
MAXCYC_SEAL ) ) {
3200 printf(
"Now: %f Target: %f \n", rhnow, rhtarget);
3203 printf(
"Cycle: %d \n", i);
3205 keepgo =
movepix(natonce, ph1id, ph2id);
3207 if ( keepgo == 1 ) {
3216 for ( j = 0; j <=
nsph; j++ ) {
3218 curvsum1 += ( j *
nsolid [ j ] );
3219 pixsum2 +=
nair [ j ];
3220 curvsum2 += ( j *
nair [ j ] );
3223 avecurv1 = ( float ) curvsum1 / (
float ) pixsum1;
3224 avecurv2 = ( float ) curvsum2 / (
float ) pixsum2;
3225 printf(
"Ave. solid curvature: %f \n", avecurv1);
3226 printf(
"Ave. air curvature: %f \n", avecurv2);
3233 int valin, ix,
iy, iz;
3234 int ix1, iy1, iz1, k;
3235 long int voltot, surftot;
3237 for ( ix = 0; ix <= 42; ix++ ) {
3242 for ( iz = 1; iz <=
SYSIZE; iz++ ) {
3243 for ( iy = 1; iy <=
SYSIZE; iy++ ) {
3244 for ( ix = 1; ix <=
SYSIZE; ix++ ) {
3245 valin =
mask [ ix ] [
iy ] [ iz ];
3252 for ( iz = 1; iz <=
SYSIZE; iz++ ) {
3253 for ( iy = 1; iy <=
SYSIZE; iy++ ) {
3254 for ( ix = 1; ix <=
SYSIZE; ix++ ) {
3255 if (
mask [ ix ] [ iy ] [ iz ] != 0 ) {
3256 valin =
mask [ ix ] [
iy ] [ iz ];
3258 for ( k = 1; k <= 6; k++ ) {
3318 if ( ( ix1 < 1 ) || ( iy1 < 1 ) || ( iz1 < 1 ) || ( ix1 >
SYSIZE ) || ( iy1 >
SYSIZE ) || ( iz1 >
SYSIZE ) ) {
3319 printf(
"%d %d %d \n", ix1, iy1, iz1);
3323 if (
mask [ ix1 ] [ iy1 ] [ iz1 ] == 0 ) {
3333 printf(
"Phase Volume Surface Volume Surface \n");
3334 printf(
" ID count count fraction fraction \n");
3341 printf(
" %d %8ld %8ld \n", k, volume [ 0 ], surface [ 0 ]);
3343 for ( k = 1; k <= 4; k++ ) {
3344 printf(
" %d %8ld %8ld %.5f %.5f\n", k, volume [ k ], surface [ k ],
3345 (
float ) volume [ k ] / (
float ) voltot, (
float ) surface [ k ] / (
float ) surftot);
3348 printf(
"Total %8ld %8ld\n\n\n", voltot, surftot);
3350 for ( k = 5; k <= 11; k++ ) {
3351 printf(
" %d %8ld %8ld\n", k, volume [ k ], surface [ k ]);
3354 printf(
" 20 %8ld %8ld\n", volume [ 20 ], surface [ 20 ]);
3356 for ( k = 24; k <= 27; k++ ) {
3357 printf(
" %d %8ld %8ld\n", k, volume [ k ], surface [ k ]);
3360 printf(
" 28 %8ld %8ld\n", volume [ 28 ], surface [ 28 ]);
3368 float s2, ss, sdiff, xtmp, ytmp;
3371 double ***normm, ***res;
3376 float *s, *xr, *sum;
3378 double t1, t2, x1, x2, u1, u2, xrad, resmax, resmin;
3379 float xtot, filval, radius, sect, sumtot, vcrit;
3380 int valin, r1, r2, i1, i2, i3, i, j, k, j1, k1;
3381 int ido, iii, jjj, ix,
iy, iz, index;
3383 char tempstr [ 256 ];
3390 s =
new float [ 61 ];
3391 xr =
new float [ 61 ];
3392 sum =
new float [ 502 ];
3399 t1 = 2. *
M_PI * u2;
3400 t2 = sqrt( -2. * log(u1) );
3403 normm [ i1 ] [ i2 ] [ i3 ] = x1;
3414 normm [ i1 ] [ i2 ] [ i3 ] = x2;
3428 F->get_next_line_in_section(0, (
long & )ido);
3436 printf(
"Number of points in correlation file is %d \n", ido);
3439 for ( i = 1; i <= ido; i++ ) {
3441 F->get_next_line_in_section(0, tempstr);
3442 sscanf(tempstr,
"%d %f", & valin, & val2);
3450 s [ i ] = ( float ) val2;
3451 xr [ i ] = ( float ) r [ i ];
3458 for ( i = 0; i < 31; i++ ) {
3460 for ( j = 0; j < 31; j++ ) {
3462 for ( k = 0; k < 31; k++ ) {
3463 xtmp = ( float ) ( iii + jjj + k * k );
3464 radius = sqrt(xtmp);
3465 r1 = ( int ) radius + 1;
3467 if ( s [ r1 ] < 0.0 ) {
3468 printf(
"%d and %d %f and %f with xtmp of %f\n", r1, r2, s [ r1 ], s [ r2 ], xtmp);
3473 xrad = radius + 1 - r1;
3474 filval = s [ r1 ] + ( s [ r2 ] - s [ r1 ] ) * xrad;
3475 filter [ i + 1 ] [ j + 1 ] [ k + 1 ] = ( filval - s2 ) / sdiff;
3484 for ( i = 1; i <=
SYSIZE; i++ ) {
3485 for ( j = 1; j <=
SYSIZE; j++ ) {
3486 for ( k = 1; k <=
SYSIZE; k++ ) {
3487 res [ i ] [ j ] [ k ] = 0.0;
3488 if ( (
float )
mask [ i ] [ j ] [ k ] == phasein ) {
3489 for ( ix = 1; ix <= 31; ix++ ) {
3499 for ( iy = 1; iy <= 31; iy++ ) {
3509 for ( iz = 1; iz <= 31; iz++ ) {
3519 res [ i ] [ j ] [ k ] += normm [ i1 ] [ j1 ] [ k1 ] * filter [ ix ] [
iy ] [ iz ];
3524 if ( res [ i ] [ j ] [ k ] > resmax ) {
3525 resmax = res [ i ] [ j ] [ k ];
3528 if ( res [ i ] [ j ] [ k ] < resmin ) {
3529 resmin = res [ i ] [ j ] [ k ];
3540 printf(
"%d out of %d\n", i,
SYSIZE);
3548 sect = ( resmax - resmin ) / 500.;
3549 for ( i = 1; i <= 500; i++ ) {
3554 for ( i = 1; i <=
SYSIZE; i++ ) {
3555 for ( j = 1; j <=
SYSIZE; j++ ) {
3556 for ( k = 1; k <=
SYSIZE; k++ ) {
3557 if ( (
float )
mask [ i ] [ j ] [ k ] == phasein ) {
3559 index = 1 + ( int ) ( ( res [ i ] [ j ] [ k ] - resmin ) / sect );
3560 if ( index > 500 ) {
3564 sum [ index ] += 1.0;
3571 sumtot = vcrit = 0.0;
3573 for ( i = 1; ( ( i <= 500 ) && ( done == 0 ) ); i++ ) {
3574 sumtot += sum [ i ] / xtot;
3575 if ( sumtot > xpt ) {
3577 vcrit = resmin + ( resmax - resmin ) * ( ytmp - 0.5 ) / 500.;
3583 printf(
"Critical volume fraction is %f\n", vcrit);
3586 for ( k = 1; k <=
SYSIZE; k++ ) {
3587 for ( j = 1; j <=
SYSIZE; j++ ) {
3588 for ( i = 1; i <=
SYSIZE; i++ ) {
3589 if ( (
float )
mask [ i ] [ j ] [ k ] == phasein ) {
3590 if ( res [ i ] [ j ] [ k ] > vcrit ) {
3591 mask [ i ] [ j ] [ k ] = phaseout;
3611 int i, j, k, alumval, alum2, valin;
3613 double volin, volf [ 5 ], surff [ 5 ], rhtest, rdesire;
3617 FILE *outfile_img = NULL, *outfile_id = NULL;
3627 printf(
"%d\n", *
seed);
3666 for ( i = 1; i <= 4; i++ ) {
3668 F->get_next_line_in_section(0, volin);
3691 printf(
"Volume %f\n", volf [ i ]);
3694 surff [ i ] = volin;
3696 printf(
"Surface %f\n", surff [ i ]);
3701 if ( ( infile = fopen(filen,
"r") ) != NULL ) {
3702 for ( k = 1; k <=
SYSIZE; k++ ) {
3703 for ( j = 1; j <=
SYSIZE; j++ ) {
3704 for ( i = 1; i <=
SYSIZE; i++ ) {
3705 if ( fscanf(infile,
"%d", & valin) != 1 ) {
3709 mask [ i ] [ j ] [ k ] = valin;
3717 for ( k = 1; k <=
SYSIZE; k++ ) {
3718 for ( j = 1; j <=
SYSIZE; j++ ) {
3719 for ( i = 1; i <=
SYSIZE; i++ ) {
3720 mask [ i ] [ j ] [ k ] =
cemreal [ i ] [ j ] [ k ];
3729 volin = volf [ 1 ] + volf [ 2 ];
3730 if ( volin < 1.0 ) {
3731 rand3d(1, alumval, volin);
3735 rdesire = ( surff [ 1 ] + surff [ 2 ] ) * (
float ) (
surface [ 1 ] +
surface [ alumval ] );
3736 if ( rdesire != 0.0 ) {
3737 if ( (
int ) rdesire <
surface [ 1 ] ) {
3738 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 1 ] ) / rdesire;
3741 rdesire = ( surff [ 3 ] + surff [ 4 ] ) * (
float ) (
surface [ 1 ] +
surface [ alumval ] );
3742 if ( rdesire != 0.0 ) {
3743 rhtest = ( 6. / 4. ) * (
float ) (
volume [ alumval ] ) / rdesire;
3751 if ( ( volf [ 1 ] + volf [ 2 ] ) > 0.0 ) {
3752 volin = volf [ 1 ] / ( volf [ 1 ] + volf [ 2 ] );
3753 if ( volin < 1.0 ) {
3758 rdesire = ( surff [ 1 ] / ( surff [ 1 ] + surff [ 2 ] ) ) * ( float ) (
surface [ 1 ] +
surface [ 2 ] );
3759 if ( rdesire != 0.0 ) {
3760 if ( (
int ) rdesire <
surface [ 1 ] ) {
3761 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 1 ] ) / rdesire;
3764 rdesire = ( surff [ 2 ] / ( surff [ 1 ] + surff [ 2 ] ) ) * ( float ) (
surface [ 1 ] +
surface [ 2 ] );
3765 if ( rdesire != 0.0 ) {
3766 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 2 ] ) / rdesire;
3775 if ( alumval == 4 ) {
3776 volin = volf [ 4 ] / ( volf [ 4 ] + volf [ 3 ] );
3779 volin = volf [ 3 ] / ( volf [ 4 ] + volf [ 3 ] );
3783 if ( volin < 1.0 ) {
3784 rand3d(alumval, alum2, volin);
3788 if ( alumval == 4 ) {
3789 rdesire = ( surff [ 4 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) (
surface [ 3 ] +
surface [ 4 ] );
3790 if ( rdesire != 0.0 ) {
3791 if ( (
int ) rdesire <
surface [ 4 ] ) {
3792 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 4 ] ) / rdesire;
3795 rdesire = ( surff [ 3 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) (
surface [ 3 ] +
surface [ 4 ] );
3796 if ( rdesire != 0.0 ) {
3797 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 3 ] ) / rdesire;
3803 rdesire = ( surff [ 3 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) (
surface [ 3 ] +
surface [ 4 ] );
3804 if ( rdesire != 0.0 ) {
3805 if ( (
int ) rdesire <
surface [ 3 ] ) {
3806 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 3 ] ) / rdesire;
3809 rdesire = ( surff [ 4 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) (
surface [ 3 ] +
surface [ 4 ] );
3810 if ( rdesire != 0.0 ) {
3811 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 4 ] ) / rdesire;
3823 F->get_value(14, (
long & )output_img);
3832 F->get_value(15, filen);
3837 if ( ( outfile_img = fopen(filen,
"w") ) == NULL ) {
3838 printf(
"Output img file %s can not be created (file %s, line %d)\n", filen, __FILE__, __LINE__);
3843 F->get_value(16, filen);
3848 if ( ( outfile_id = fopen(filen,
"w") ) == NULL ) {
3849 printf(
"Output id file %s can not be created (file %s, line %d)\n", filen, __FILE__, __LINE__);
3854 for ( k = 0; k <
SYSIZE; k++ ) {
3855 for ( j = 0; j <
SYSIZE; j++ ) {
3856 for ( i = 0; i <
SYSIZE; i++ ) {
3857 mic [ i ] [ j ] [ k ] =
mask [ i + 1 ] [ j + 1 ] [ k + 1 ];
3858 micorig [ i ] [ j ] [ k ] =
mic [ i ] [ j ] [ k ];
3860 fprintf(outfile_img,
"%d\n",
mic [ i ] [ j ] [ k ]);
3861 fprintf(outfile_id,
"%ld\n",
micpart [ i ] [ j ] [ k ]);
3868 fclose(outfile_img);
3888 double slagin, CHperslag;
3891 for ( i = 0; i <=
EMPTYP; i++ ) {
4225 if ( ( slagfile = fopen(
"slagchar.dat",
"r") ) == NULL ) {
4230 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4234 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4238 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4243 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4248 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4253 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4258 if ( fscanf(slagfile,
"%lf", &
slagcasi) != 1 ) {
4262 if ( fscanf(slagfile,
"%lf", &
slaghydcasi) != 1 ) {
4266 if ( fscanf(slagfile,
"%lf", &
siperslag) != 1 ) {
4270 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4275 if ( fscanf(slagfile,
"%lf", &
slagc3a) != 1 ) {
4279 if ( fscanf(slagfile,
"%lf", &
slagreact) != 1 ) {
4303 if ( CHperslag < 0.0 ) {
4314 printf(
"Error in range of C3A/slag value... reset to 1.0 \n");
4324 int edgeback, x2, y2, z2;
4331 for ( ip = 0; ( ( ip <
NEIGHBORS ) && ( edgeback == 0 ) ); ip++ ) {
4332 x2 = xck +
xoff [ ip ];
4333 y2 = yck +
yoff [ ip ];
4334 z2 = zck +
zoff [ ip ];
4364 return ( edgeback );
4373 int i, xid, yid, zid, phid, edgef, phread, cshcyc;
4381 for ( i = low; i <= high; i++ ) {
4386 for ( xid = 0; xid <
SYSIZE; xid++ ) {
4387 for ( yid = 0; yid <
SYSIZE; yid++ ) {
4388 for ( zid = 0; zid <
SYSIZE; zid++ ) {
4389 phread =
mic [ xid ] [ yid ] [ zid ];
4391 if ( ( cshexflag == 1 ) && ( phread ==
CSH ) ) {
4392 cshcyc =
cshage [ xid ] [ yid ] [ zid ];
4399 for ( i = low; ( ( i <= high ) && ( phid == 60 ) ); i++ ) {
4400 if (
mic [ xid ] [ yid ] [ zid ] == i ) {
4415 else if ( i ==
C3S ) {
4417 }
else if ( i ==
C2S ) {
4419 }
else if ( i ==
C3A ) {
4421 }
else if ( i ==
C4AF ) {
4423 }
else if ( i ==
GYPSUM ) {
4431 }
else if ( i ==
POZZ ) {
4433 }
else if ( i ==
SLAG ) {
4435 }
else if ( i ==
ETTR ) {
4446 if ( ( cycid != 0 ) && (
soluble [ phid ] == 1 ) ) {
4471 int effort, tries, xmod, ymod, zmod;
4472 struct ants *antnew;
4477 while ( ( effort == 0 ) && ( tries < 500 ) ) {
4479 xmod = ( -extent ) + (
int ) ( ( 2. * ( float ) extent + 1. ) *
ran1(
seed) );
4480 ymod = ( -extent ) + (
int ) ( ( 2. * ( float ) extent + 1. ) *
ran1(
seed) );
4481 zmod = ( -extent ) + (
int ) ( ( 2. * ( float ) extent + 1. ) *
ran1(
seed) );
4482 if ( xmod > extent ) {
4486 if ( ymod > extent ) {
4490 if ( zmod > extent ) {
4500 }
else if ( xmod < 0 ) {
4506 }
else if ( zmod < 0 ) {
4512 }
else if ( ymod >=
SYSIZE ) {
4516 if (
mic [ xmod ] [ ymod ] [ zmod ] ==
POROSITY ) {
4522 antnew = (
struct ants * ) malloc(
sizeof(
struct ants ) );
4545 int nfound, ix,
iy, iz, qxlo, qxhi, qylo, qyhi, qzlo, qzhi;
4546 int hx, hy, hz, boxhalf;
4548 boxhalf = boxsize / 2;
4550 qxlo = qx - boxhalf;
4551 qxhi = qx + boxhalf;
4552 qylo = qy - boxhalf;
4553 qyhi = qy + boxhalf;
4554 qzlo = qz - boxhalf;
4555 qzhi = qz + boxhalf;
4558 for ( ix = qxlo; ix <= qxhi; ix++ ) {
4562 }
else if ( hx >=
SYSIZE ) {
4566 for ( iy = qylo; iy <= qyhi; iy++ ) {
4570 }
else if ( hy >=
SYSIZE ) {
4574 for ( iz = qzlo; iz <= qzhi; iz++ ) {
4578 }
else if ( hz >=
SYSIZE ) {
4583 if ( (
mic [ hx ] [ hy ] [ hz ] <
C3S ) || (
mic [ hx ] [ hy ] [ hz ] >
ABSGYP ) ) {
4601 int nfound, ix,
iy, iz, qxlo, qxhi, qylo, qyhi, qzlo, qzhi;
4602 int hx, hy, hz, boxhalf;
4604 boxhalf = boxsize / 2;
4606 qxlo = qx - boxhalf;
4607 qxhi = qx + boxhalf;
4608 qylo = qy - boxhalf;
4609 qyhi = qy + boxhalf;
4610 qzlo = qz - boxhalf;
4611 qzhi = qz + boxhalf;
4614 for ( ix = qxlo; ix <= qxhi; ix++ ) {
4618 }
else if ( hx >=
SYSIZE ) {
4622 for ( iy = qylo; iy <= qyhi; iy++ ) {
4626 }
else if ( hy >=
SYSIZE ) {
4630 for ( iz = qzlo; iz <= qzhi; iz++ ) {
4634 }
else if ( hz >=
SYSIZE ) {
4639 if ( (
mic [ hx ] [ hy ] [ hz ] <
C3S ) || (
mic [ hx ] [ hy ] [ hz ] >
POZZ ) ) {
4658 struct togo *headtogo, *tailtogo, *newtogo, *lasttogo, *onetogo;
4660 int px, py, pz, placed, cntpore, cntmax;
4663 headtogo = (
struct togo * ) malloc(
sizeof(
struct togo ) );
4664 headtogo->
x = headtogo->
y = headtogo->
z = ( -1 );
4665 headtogo->
npore = 0;
4668 tailtogo = headtogo;
4672 printf(
"In makeinert with %ld needed elements \n", ndesire);
4677 for ( idesire = 2; idesire <= ndesire; idesire++ ) {
4678 newtogo = (
struct togo * ) malloc(
sizeof(
struct togo ) );
4680 newtogo->
x = newtogo->
y = newtogo->
z = ( -1 );
4687 for ( px = 0; px <
SYSIZE; px++ ) {
4688 for ( py = 0; py <
SYSIZE; py++ ) {
4689 for ( pz = 0; pz <
SYSIZE; pz++ ) {
4692 if ( cntpore > cntmax ) {
4698 if ( cntpore > ( tailtogo->
npore ) ) {
4700 lasttogo = tailtogo;
4701 while ( placed == 0 ) {
4703 if ( newtogo == NULL ) {
4706 if ( cntpore <= ( newtogo->
npore ) ) {
4711 if ( placed == 0 ) {
4716 onetogo = (
struct togo * ) malloc(
sizeof(
struct togo ) );
4720 onetogo->
npore = cntpore;
4722 if ( placed == 2 ) {
4729 if ( placed == 1 ) {
4737 lasttogo = tailtogo;
4751 for ( idesire = 1; idesire <= ndesire; idesire++ ) {
4755 if ( px != ( -1 ) ) {
4761 lasttogo = headtogo;
4782 int check, sump, xchr, ychr, zchr, fchr, i1, action, numnear;
4784 int mstest = 0, mstest2 = 0;
4793 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
4799 sump *=
moveone(& xchr, & ychr, & zchr, & action, sump);
4800 if ( action == 0 ) {
4801 printf(
"Error in value of action in extpozz \n");
4804 check =
mic [ xchr ] [ ychr ] [ zchr ];
4807 if ( xchr != xpres ) {
4812 if ( ychr != ypres ) {
4817 if ( zchr != zpres ) {
4824 if ( (
faces [ xpres ] [ ypres ] [ zpres ] == 0 ) || ( mstest ==
faces [ xpres ] [ ypres ] [ zpres ] ) || ( mstest2 ==
faces [ xpres ] [ ypres ] [ zpres ] ) ) {
4826 faces [ xchr ] [ ychr ] [ zchr ] =
faces [ xpres ] [ ypres ] [ zpres ];
4837 while ( fchr == 0 ) {
4855 check =
mic [ xchr ] [ ychr ] [ zchr ];
4861 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
4876 int nc3aext, ncshext, nchext, ngypext, nanhext, plok;
4877 int nsum5, nsum4, nsum3, nsum2, nhemext, nsum6, nc4aext;
4878 int phid, phnew, plnew, cread;
4879 int i, xloop, yloop, zloop, xc, yc;
4883 int placed, cshrand, maxsulfate, msface;
4884 long int ncshgo, nsurf, suminit;
4885 long int xext, nhgd, npchext, nslagc3a = 0;
4886 float pdis, plfh3, fchext, fc3aext, fanhext, mass_now, mass_fa_now, tot_mass, heatfill;
4887 float dfact, dfact1, molesdh2o, h2oinit, heat4, fhemext, fc4aext;
4888 float pconvert, pc3scsh, pc2scsh, calcx, calcy, calcz, tdisfact;
4889 float frafm, frettr, frhyg, frtot, mc3ar, mc4ar, p3init;
4890 struct ants *antadd;
4894 npchext = ncshgo = cshrand = 0;
4900 for ( i = 0; i <=
EMPTYP; i++ ) {
4911 printf(
"Returned from passone \n");
4988 heatfill = ( float ) ( count [
INERT ] + count [
SLAG ] + count [ POZZ ] + count [
CACL2 ] + count [
ASG ] + count [
CAS2 ] ) /
cemmass;
5006 printf(
"Calculated w/c is %.4f\n",
w_to_c);
5007 printf(
"Calculated s/c is %.4f \n",
s_to_c);
5008 printf(
"Calculated heat conversion factor is %f \n",
heat_cf);
5009 printf(
"Calculated mass fractions of water and filler are %.4f and %.4f \n",
5026 printf(
"ctest is %ld\n", ctest);
5029 if ( (
float ) ctest > ( 2.5 * (
float ) ( count [
DIFFC3A ] + count [
DIFFC4A ] ) ) ) {
5030 ctest = (
long int ) ( 2.5 * (
float ) ( count [
DIFFC3A ] + count [
DIFFC4A ] ) );
5033 for ( i = 0; i <=
EMPTYP; i++ ) {
5046 if ( i == DIFFC4A ) {
5057 }
else if ( i ==
DIFFCH ) {
5073 }
else if ( i ==
DIFFAS ) {
5091 }
else if ( i ==
C3S ) {
5093 mass_now +=
specgrav [
C3S ] * ( float ) count [ C3S ];
5095 }
else if ( i ==
C2S ) {
5097 mass_now +=
specgrav [
C2S ] * ( float ) count [ C2S ];
5099 }
else if ( i ==
C3A ) {
5101 mass_now +=
specgrav [
C3A ] * ( float ) count [ C3A ];
5103 mc4ar = ( c4afinit - ( float ) count [
C4AF ] ) /
molarv [
C4AF ];
5104 if ( ( mc3ar + mc4ar ) > 0.0 ) {
5105 frhyg = ( mc3ar / ( mc3ar + mc4ar ) ) * ( float ) count [
C3AH6 ] /
molarv [
C3AH6 ];
5112 frtot = frafm + frettr + frhyg;
5113 if ( frtot > 0.0 ) {
5121 }
else if ( i ==
C4AF ) {
5122 alpha += ( float ) ( c4afinit - count [
C4AF ] );
5123 mass_now +=
specgrav [
C4AF ] * ( float ) count [ C4AF ];
5125 mc4ar = ( c4afinit - ( float ) count [ C4AF ] ) /
molarv [
C4AF ];
5126 if ( ( mc3ar + mc4ar ) > 0.0 ) {
5127 frhyg = ( mc4ar / ( mc3ar + mc4ar ) ) * ( float ) count [
C3AH6 ] /
molarv [
C3AH6 ];
5133 frtot = frettr + frhyg;
5134 if ( frtot > 0.0 ) {
5137 heat4 += frhyg * .418 * ( float ) ( c4afinit - count [ C4AF ] ) *
specgrav [
C4AF ];
5138 heat4 += frettr * .725 * ( float ) ( c4afinit - count [ C4AF ] ) *
specgrav [
C4AF ];
5166 if ( suminit != 0 ) {
5210 fprintf(
heatfile,
"Cycle time(h) alpha_vol alpha_mass heat4(kJ/kg_solid) Gsratio2 G-s_ratio\n");
5225 printf(
"All water consumed at cycle %d \n",
cyccnt);
5248 fprintf(
phasfile,
"#Cycle DoH Time Porosity C3S C2S C3A C4AF GYPSUM HEMIHYD ANHYDRITE POZZ INERT SLAG ASG CAS2 CH CSH C3AH6 ETTR ETTRC4AF AFM FH3 POZZCSH SLAGCSH CACL2 FREIDEL STRAT GYPSUMS CACO3 AFMC AGG ABSGYP EMPTYP HDCSH water_left \n");
5249 fprintf(
perc_phases,
"#Cyc Time[h] DoH SOL_per SOL_tot| Porosity C3S C2S C3A C4AF GYPSUM HEMIHYD ANHYDRITE POZZ INERT SLAG ASG CAS2 CH CSH C3AH6 ETTR ETTRC4AF AFM FH3 POZZCSH SLAGCSH CACL2 FREIDEL STRAT GYPSUMS CACO3 AFMC AGG ABSGYP EMPTYP HDCSH\n");
5259 for ( i = 0; i <=
HDCSH; i++ ) {
5264 fprintf(
phasfile,
"%ld ", count [ i ]);
5308 ( ( ( (
float ) count [
GYPSUM ] + 1.42 * (
float ) count [
ANHYDRITE ] + 1.4 *
5310 1.42 * (
float )
anhinit + 1.4 * (
float )
heminit + ( (
float )
netbar / 3.30 ) ) ) < 0.25 ) ) ) {
5313 printf(
"Ettringite is soluble beginning at cycle %d \n", cycle);
5364 printf(
"CH dissolution probability goes from %f ",
disprob [
CH ]);
5372 printf(
"to %f \n",
disprob [ CH ]);
5384 1.42 * (
float )
anhinit + 1.4 * (
float )
heminit ) * 0.05 ) ) ||
5391 maxsulfate = count [
DIFFGYP ];
5402 if ( maxsulfate > 0 ) {
5415 if ( (
soluble [
C3S ] == 0 ) && ( ( cycle > 1 ) || ( count [
ETTR ] > 0 ) || ( count [
AFM ] > 0 ) || ( count [
ETTRC4AF ] > 0 ) ) ) {
5497 if (
disprob [ C3S ] == disbase [ C3S ] ) {
5506 if (
disprob [ C3S ] == disbase [ C3S ] ) {
5515 if (
disprob [ C3S ] == disbase [ C3S ] ) {
5526 printf(
"Silicate probabilities: %f %f\n",
disprob [ C3S ],
disprob [ C2S ]);
5557 if (
disprob [ GYPSUMS ] > ( disbase [ GYPSUMS ] ) ) {
5727 for ( xloop = 0; xloop <
SYSIZE; xloop++ ) {
5728 for ( yloop = 0; yloop <
SYSIZE; yloop++ ) {
5729 for ( zloop = 0; zloop <
SYSIZE; zloop++ ) {
5730 if (
mic [ xloop ] [ yloop ] [ zloop ] >
OFFSET ) {
5731 phid =
mic [ xloop ] [ yloop ] [ zloop ] -
OFFSET;
5734 if ( ( plnew < 0 ) || ( plnew >=
NEIGHBORS ) ) {
5738 xc = xloop +
xoff [ plnew ];
5739 yc = yloop +
yoff [ plnew ];
5740 zc = zloop +
zoff [ plnew ];
5749 if ( xc >= SYSIZE ) {
5753 if ( yc >= SYSIZE ) {
5761 if ( zc >= SYSIZE ) {
5773 count [ phid ] -= 1;
5775 if ( phid ==
C3AH6 ) {
5780 if ( phid == C4AF ) {
5782 if ( ( plfh3 < 0.0 ) || ( plfh3 > 1.0 ) ) {
5788 if ( plfh3 <= 0.5453 ) {
5801 count [ phnew ] += 1;
5802 mic [ xc ] [ yc ] [ zc ] = phnew;
5803 antadd = (
struct ants * ) malloc(
sizeof(
struct ants ) );
5817 if ( ( phid == C3S ) || ( phid ==
C2S ) ) {
5819 if ( ( ( phid == C2S ) && ( plfh3 <= pc2scsh ) ) || ( plfh3 <= pc3scsh ) ) {
5831 if ( placed != 0 ) {
5840 if ( ( phid == C2S ) && ( pc2scsh > 1.0 ) ) {
5842 if ( plfh3 <= ( pc2scsh - 1.0 ) ) {
5854 if ( placed != 0 ) {
5863 mic [ xloop ] [ yloop ] [ zloop ] -=
OFFSET;
5872 if ( ( count [
POZZ ] >= 13000 ) && (
chnew < ( 0.15 * SYSIZE * SYSIZE * SYSIZE ) ) && (
csh2flag == 1 ) ) {
5873 if (
mic [ xloop ] [ yloop ] [ zloop ] ==
CSH ) {
5874 if ( (
countbox(3, xloop, yloop, zloop) ) >= 1 ) {
5884 cycnew =
cshage [ xloop ] [ yloop ] [ zloop ];
5886 if ( calcy > 1.0 ) {
5887 calcz = calcy - 1.0;
5889 printf(
"Problem of not creating enough pozzolanic CSH during CSH conversion \n");
5890 printf(
"Current temperature is %f C\n",
temp_cur);
5893 if ( plfh3 <= calcy ) {
5897 mic [ xloop ] [ yloop ] [ zloop ] =
DIFFCH;
5902 antadd = (
struct ants * ) malloc(
sizeof(
struct ants ) );
5927 calcx = ( 19.86 /
molarvcsh [ cycnew ] ) - ( 1. - calcy );
5929 if ( plfh3 < calcx ) {
5938 if (
mic [ xloop ] [ yloop ] [ zloop ] == SLAG ) {
5939 if ( (
countbox(3, xloop, yloop, zloop) ) >= 1 ) {
5943 count [
SLAG ] -= 1;
5956 msface = ( int ) ( 3. *
ran1(
seed) + 1. );
5961 faces [ xloop ] [ yloop ] [ zloop ] = msface;
5967 mic [ xloop ] [ yloop ] [ zloop ] =
EMPTYP;
5977 while ( p3init > 1.0 ) {
5983 if ( plfh3 < p3init ) {
5999 if ( ncshgo != 0 ) {
6000 printf(
"CSH dissolved is %ld \n", ncshgo);
6003 if ( npchext > 0 ) {
6004 printf(
"npchext is %ld at cycle %d \n", npchext, cycle);
6011 if ( cshrand != 0 ) {
6012 printf(
"cshrand is %d \n", cshrand);
6026 nchext = ( int ) fchext;
6027 if ( fchext > (
float ) nchext ) {
6029 if ( ( fchext - (
float ) nchext ) > pdis ) {
6048 nc3aext = ( int ) ( fc3aext + nslagc3a );
6049 if ( fc3aext > (
float ) nc3aext ) {
6051 if ( ( fc3aext - (
float ) nc3aext ) > pdis ) {
6056 fc4aext = 0.696 * ( float )
discount [ C4AF ];
6057 nc4aext = ( int ) fc4aext;
6058 if ( fc4aext > (
float ) nc4aext ) {
6060 if ( ( fc4aext - (
float ) nc4aext ) > pdis ) {
6072 fanhext = ( float ) discount [
ANHYDRITE ];
6073 nanhext = ( int ) fanhext;
6074 if ( fanhext > (
float ) nanhext ) {
6076 if ( ( fanhext - (
float ) nanhext ) > pdis ) {
6085 fhemext = ( float ) discount [
HEMIHYD ];
6088 nhemext = ( int ) fhemext;
6089 if ( fhemext > (
float ) nhemext ) {
6091 if ( ( fhemext - (
float ) nhemext ) > pdis ) {
6099 count [
DIFFCH ] += nchext;
6104 nsum2 = nchext + ncshext;
6105 nsum3 = nsum2 + nc3aext;
6106 nsum4 = nsum3 + nc4aext;
6107 nsum5 = nsum4 + ngypext;
6108 nsum6 = nsum5 + nhemext;
6112 for ( xext = 1; xext <= ( nsum6 + nanhext ); xext++ ) {
6116 xc = ( int ) ( (
float ) SYSIZE *
ran1(
seed) );
6117 yc = ( int ) ( (
float ) SYSIZE *
ran1(
seed) );
6118 zc = ( int ) ( (
float ) SYSIZE *
ran1(
seed) );
6119 if ( xc >= SYSIZE ) {
6123 if ( yc >= SYSIZE ) {
6127 if ( zc >= SYSIZE ) {
6137 if ( xext > nsum6 ) {
6139 }
else if ( xext > nsum5 ) {
6141 }
else if ( xext > nsum4 ) {
6143 }
else if ( xext > nsum3 ) {
6145 }
else if ( xext > nsum2 ) {
6147 }
else if ( xext > nchext ) {
6151 mic [ xc ] [ yc ] [ zc ] = phid;
6154 antadd = (
struct ants * ) malloc(
sizeof(
struct ants ) );
6170 printf(
"Could not place extra diffusing species (too few POROSITY voxels), continuing (line %d),\n", __LINE__);
6174 }
while ( plok == 0 );
6180 printf(
"Dissolved- %ld %ld %ld %ld %ld %ld %ld %ld %ld %ld %ld %ld\n", count [
DIFFCSH ],
6203 printf(
"C3AH6 dissolved- %ld with prob. of %f \n", nhgd,
disprob [
C3AH6 ]);
6220 int success, cpores;
6223 for ( ic = 1; ic <= nneed; ic++ ) {
6225 while ( success == 0 ) {
6242 if ( ( randid !=
CACO3 ) && ( randid !=
INERT ) ) {
6243 mic [ ix ] [
iy ] [ iz ] = randid;
6248 if ( cpores >= 26 ) {
6249 mic [ ix ] [
iy ] [ iz ] = randid;
6262 int sx, sy, sz, jx = 0, jy = 0, jz = 0, faceid;
6264 for ( sx = 0; sx <
SYSIZE; sx++ ) {
6265 for ( sy = 0; sy <
SYSIZE; sy++ ) {
6266 for ( sz = 0; sz <
SYSIZE; sz++ ) {
6268 for ( faceid = 0; faceid < 6; faceid++ ) {
6269 if ( faceid == 1 ) {
6277 }
else if ( faceid == 0 ) {
6279 if ( jx > ( SYSIZE - 1 ) ) {
6285 }
else if ( faceid == 2 ) {
6287 if ( jy > ( SYSIZE - 1 ) ) {
6293 }
else if ( faceid == 3 ) {
6301 }
else if ( faceid == 4 ) {
6303 if ( jz > ( SYSIZE - 1 ) ) {
6309 }
else if ( faceid == 5 ) {
6320 if ( (
mic [ jx ] [ jy ] [ jz ] ==
C3S ) || (
mic [ jx ] [ jy ] [ jz ] ==
C2S ) || (
mic [ jx ] [ jy ] [ jz ] ==
C3A ) || (
mic [ jx ] [ jy ] [ jz ] ==
C4AF ) || (
mic [ jx ] [ jy ] [ jz ] ==
INERT ) || (
mic [ jx ] [ jy ] [ jz ] ==
CACO3 ) ) {
6322 if ( (
mic [ jx ] [ jy ] [ jz ] ==
C3S ) || (
mic [ jx ] [ jy ] [ jz ] ==
C2S ) || (
mic [ jx ] [ jy ] [ jz ] ==
C3A ) || (
mic [ jx ] [ jy ] [ jz ] ==
C4AF ) ) {
6333 printf(
"Cement surface count is %ld \n",
scntcement);
6334 printf(
"Total surface count is %ld \n",
scnttotal);
6340 printf(
"Surface fraction is %f \n",
surffract);
6351 long int nresat = 0;
6353 for ( sx = 0; sx <
SYSIZE; sx++ ) {
6354 for ( sy = 0; sy <
SYSIZE; sy++ ) {
6355 for ( sz = 0; sz <
SYSIZE; sz++ ) {
6356 if (
mic [ sx ] [ sy ] [ sz ] ==
EMPTYP ) {
6368 printf(
"Number resaturated is %ld \n", nresat);
6377 char extension [ 10 ];
6380 prefix = (
char * ) malloc(80);
6388 sprintf(extension,
"%04d",
icyc);
6389 strcpy(prefix,
"unperc/out5.");
6390 strcat(prefix, extension);
6391 strcat(prefix,
".img");
6393 printf(
"Name of output file is %s\n", prefix);
6396 if ( ( imgout = fopen(prefix,
"w") ) == NULL ) {
6397 printf(
"File %s can not be opened\n", prefix);
6402 for (
int dz = 0; dz <
SYSIZE; dz++ ) {
6403 for (
int dy = 0; dy <
SYSIZE; dy++ ) {
6404 for (
int dx = 0; dx <
SYSIZE; dx++ ) {
6405 fprintf(imgout,
"%d\n", m [ dx ] [ dy ] [ dz ]);
6412 printf(
"Unpercolated file %s wrote\n", prefix);
6420 int fidc3s, fidc2s, fidc3a, fidc4af, fidgyp, fidagg, ffac3a;
6421 int fidhem, fidanh, fidcaco3;
6424 int ix,
iy, iz, phtodo;
6436 printf(
"Seed %d\n", *
seed);
6441 F->get_value(13, (
long & )read_micr);
6450 printf(
"Enter name of file to read initial microstructure from \n");
6453 F->get_value(1, filei);
6461 printf(
"%s\n", filei);
6470 printf(
"Enter IDs in file for C3S, C2S, C3A, C4AF, Gypsum, Hemihydrate, Anhydrite, Aggregate, CaCO3\n");
6485 printf(
"%d %d %d %d %d %d %d %d %d\n", fidc3s, fidc2s, fidc3a, fidc4af, fidgyp, fidhem, fidanh, fidagg, fidcaco3);
6486 printf(
"Enter ID in file for C3A in fly ash (default=35)\n");
6491 printf(
"%d\n", ffac3a);
6496 if ( ( infile = fopen(filei,
"r") ) == NULL ) {
6497 printf(
"CEMHYD3D microstructure file %s not found\n", filei);
6501 for ( iz = 0; iz <
SYSIZE; iz++ ) {
6502 for ( iy = 0; iy <
SYSIZE; iy++ ) {
6503 for ( ix = 0; ix <
SYSIZE; ix++ ) {
6505 faces [ ix ] [
iy ] [ iz ] = 0;
6507 if ( fscanf(infile,
"%ld", & valin) == EOF ) {
6508 printf(
"End of file %s reached, terminating\n", filei);
6513 printf(
"Error in the reading at x=%d y=%d z=%d, value %ld\n", ix, iy, iz, valin);
6517 mic [ ix ] [
iy ] [ iz ] = valin;
6518 if ( valin == fidc3s ) {
6520 }
else if ( valin == fidc2s ) {
6522 }
else if ( ( valin == fidc3a ) || ( valin == ffac3a ) ) {
6524 }
else if ( valin == fidc4af ) {
6526 }
else if ( valin == fidgyp ) {
6528 }
else if ( valin == fidanh ) {
6530 }
else if ( valin == fidhem ) {
6532 }
else if ( valin == fidcaco3 ) {
6534 }
else if ( valin == fidagg ) {
6546 printf(
"Generated microstructure from memory will be used\n");
6555 printf(
"Enter name of file to read particle IDs from \n");
6558 F->get_value(2, filei);
6566 printf(
"%s\n", filei);
6569 if ( ( infile = fopen(filei,
"r") ) == NULL ) {
6570 printf(
"CEMHYD3D microstructure ID file %s not found\n", filei);
6574 for ( iz = 0; iz <
SYSIZE; iz++ ) {
6575 for ( iy = 0; iy <
SYSIZE; iy++ ) {
6576 for ( ix = 0; ix <
SYSIZE; ix++ ) {
6577 if ( fscanf(infile,
"%ld", & valin) == EOF ) {
6578 printf(
"End of file %s reached, terminating\n", filei);
6590 printf(
"Generated microstructure ID from memory will be used\n");
6601 printf(
"Enter number of one pixel particles to add (0 to quit) \n");
6607 printf(
"%ld\n", nadd);
6610 while ( nadd > 0 ) {
6612 printf(
"Enter phase to add \n");
6613 printf(
" C3S 1 \n");
6614 printf(
" C2S 2 \n");
6615 printf(
" C3A 3 \n");
6616 printf(
" C4AF 4 \n");
6617 printf(
" GYPSUM 5 \n");
6618 printf(
" HEMIHYD 6 \n");
6619 printf(
" ANHYDRITE 7 \n");
6620 printf(
" POZZ 8 \n");
6621 printf(
" INERT 9 \n");
6622 printf(
" SLAG 10 \n");
6623 printf(
" ASG 11 \n");
6624 printf(
" CAS2 12 \n");
6625 printf(
" CH 13 \n");
6626 printf(
" CSH 14 \n");
6627 printf(
" C3AH6 15 \n");
6628 printf(
" Ettringite 16 \n");
6629 printf(
" Stable Ettringite from C4AF 17 \n");
6630 printf(
" AFM 18 \n");
6631 printf(
" FH3 19 \n");
6632 printf(
" POZZCSH 20 \n");
6633 printf(
" SLAGCSH 21 \n");
6634 printf(
" CACL2 22 \n");
6635 printf(
" Friedels salt 23 \n");
6636 printf(
" Stratlingite 24 \n");
6637 printf(
" Calcium carbonate 26 \n");
6642 printf(
"%d \n", phtodo);
6644 if ( ( phtodo < 0 ) || ( phtodo >
CACO3 ) ) {
6645 printf(
"Error in phase input for one pixel particles \n");
6652 printf(
"Enter number of one pixel particles to add (0 to quit) \n");
6658 printf(
"%ld\n", nadd);
6666 F->get_value(3, (
long & )
sealed);
6676 F->get_value(23, (
long & )
ntimes);
6686 F->get_value(24,
pnucch);
6712 F->get_value(28,
pnuchg);
6738 F->get_value(17, (
long & )
burnfreq);
6746 #ifdef __TM_MODULE //OOFEM transport module 6751 printf(
"Enter cycle frequency for checking percolation of solids (set) \n");
6752 printf(
"%d\n", burnfreq);
6757 F->get_value(18, (
long & )
setfreq);
6764 printf(
"Enter cycle frequency for checking hydration of particles \n");
6765 printf(
"%d\n", setfreq);
6788 printf(
"Enter apparent activation energy for hydration in kJ/mole \n");
6791 F->get_value(5,
E_act);
6798 printf(
"%f \n",
E_act);
6799 printf(
"Enter apparent activation energy for pozzolanic reactions in kJ/mole \n");
6810 printf(
"Enter apparent activation energy for slag reactions in kJ/mole \n");
6821 printf(
"Enter kinetic factor to convert cycles to time for 25 C \n");
6824 F->get_value(8,
beta);
6831 printf(
"%f \n",
beta);
6832 printf(
"Enter mass fraction of aggregate in concrete \n");
6843 printf(
"Enter kg of cement (+admixtures) in 1 m3\n");
6854 printf(
"Enter heat capacity of aggregate in concrete J/g/C\n");
6857 F->get_value(11,
Cp_agg);
6867 printf(
"Enter heat capacity of cement J/g/C\n");
6880 printf(
"CSH to pozzolanic CSH 0) prohibited or 1) allowed \n");
6886 printf(
"CH precipitation on aggregate surfaces 0) prohibited or 1) allowed \n");
6892 printf(
"Enter number of cycles before executing total resaturation \n");
6898 printf(
"Enter choice for C-S-H geometry 0) random or 1) plates \n");
6904 printf(
"Does pH influence hydration kinetics 0) no or 1) yes \n");
6915 F->get_value(37,
Vol_FA);
6916 F->get_value(38,
Vol_CA);
7060 printf(
"After init routine \n");
7067 fileperc = fopen(
"percpore.out",
"w");
7069 fprintf(
fileperc,
"#cycle alpha #through #tot_phase\n");
7078 strcpy(
pHname,
"pHfile.out");
7080 fprintf(
pHfile,
"Cycle time(h) alpha_mass pH sigma [Na+] [K+] [Ca++] [SO4--] activityCa activityOH activitySO4 activityK molesSyngenite\n");
7081 strcpy(fileo,
"out5.img");
7084 strcpy(
ppsname,
"percpore.out");
7085 strcpy(
phrname,
"parthydr.out");
7088 fprintf(
adiafile,
"Cyc Time(h) Temperature Alpha Krate Cp_now Mass_cem kpozz/khyd kslag/khyd DoR_Blend\n");
7090 fprintf(
elasfile,
"Cyc\tTime[h] Alpha E_paste nu_paste E_paste_fil nu_paste_fil E_mortar nu_mortar E_concrete nu_concrete\n");
7091 CSHfile = fopen(
"CSH.out",
"w");
7092 fprintf(
CSHfile,
"cyc alpha [h] HD CSH-total HD/CSH\n");
7093 system(
"mkdir perc 2> /dev/null");
7094 system(
"mkdir unperc 2> /dev/null");
7098 infoperc = fopen(
"perc/info.dat",
"w");
7099 fprintf(
infoperc,
"#Cyc DoH Time[h]\n");
7107 krate = exp( -( 1000. *
E_act / 8.314 ) * ( ( 1. / (
temp_cur + 273.15 ) ) - ( 1. / 298.15 ) ) );
7151 while ( flag == 0 ? (
time_cur < hydrationTime ) : counter < flag ) {
7217 krate = exp( -( 1000. *
E_act / 8.314 ) * ( ( 1. / (
temp_cur + 273.15 ) ) - ( 1. / 298.15 ) ) );
7257 fprintf(
adiafile,
"%d \t %.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f\n",
cyccnt,
time_cur,
temp_cur,
alpha_cur,
krate,
Cp_now,
mass_cem_now,
kpozz /
krate,
kslag /
krate,
alpha_fa_cur);
7258 fprintf(
heatfile,
"%d %f %f %f %f %f %f \n",
7284 printf(
"Switching to self-desiccating at cycle %d \n",
cyccnt);
7324 double E_paste_SCM_filler_air, nu_paste_SCM_filler_air, E_mortar, nu_mortar;
7331 fprintf(
elasfile,
"%d \t%.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f\n",
cyccnt,
time_cur,
alpha_cur,
last_values [ 2 ],
last_values [ 3 ], E_paste_SCM_filler_air, nu_paste_SCM_filler_air, E_mortar, nu_mortar,
last_values [ 4 ],
last_values [ 5 ]);
7374 #ifdef __TM_MODULE //OOFEM transport module 7377 double castingTime = 0.;
7381 if ( TargTime - castingTime <= 0 ) {
7390 printf(
"Cannot go backwards in hydration, TargTime = %f s < LastCallTime = %f s\n", TargTime,
LastCallTime);
7394 #ifdef __TM_MODULE //OOFEM transport module 7396 if ( !cemhydmat->
nowarnings.
at(4) && ( GiveTemp > 200. ) ) {
7397 printf(
"Temperature exceeds 200 C (file %s, line %d),\n", __FILE__, __LINE__);
7401 if ( GiveTemp > 200. ) {
7402 printf(
"Temperature exceeds 200 C (file %s, line %d),\n", __FILE__, __LINE__);
7411 disrealnew(GiveTemp, ( TargTime - castingTime ) / 3600., 0);
7415 printf(
"time_cur %f [h]\n",
time_cur);
7451 printf(
"PartHeat %f [W/m3 of concrete]\n",
PartHeat);
7463 double castingTime = 0.;
7503 if ( cycle > maxcyc ) {
7505 printf(
"Target DoH %f (now is %f) was not reached after %d cycles (%s, line %d)\n", DesiredDoH,
alpha_cur, maxcyc, __FILE__, __LINE__);
7549 double capacityPaste, capacityConcrete;
7551 capacityPaste *= ( 1. - 0.26 * ( 1 - pow(2.71828, -2.9 *
alpha_cur) ) );
7554 return capacityConcrete;
7571 double castingTime = 0;
7608 long int ntop, nthrough, ncur, nnew, ntot, nphc;
7610 int xl, xh, j1, k1, px, py, pz, qx, qy, qz, xcn, ycn, zcn;
7611 int x1, y1, z1, igood;
7612 int *nmatx, *nmaty, *nmatz, *nnewx, *nnewy, *nnewz;
7615 float mass_burn = 0.0, alpha_burn = 0.0;
7617 nmatx =
new int [
SIZE2D ];
7618 nmaty =
new int [
SIZE2D ];
7619 nmatz =
new int [
SIZE2D ];
7620 nnewx =
new int [
SIZE2D ];
7621 nnewy =
new int [
SIZE2D ];
7622 nnewz =
new int [
SIZE2D ];
7638 for ( k = 0; k <
SYSIZE; k++ ) {
7639 for ( j = 0; j <
SYSIZE; j++ ) {
7644 px =
cx(i, j, k, d1, d2, d3);
7645 py =
cy(i, j, k, d1, d2, d3);
7646 pz =
cz(i, j, k, d1, d2, d3);
7647 if (
mic [ px ] [ py ] [ pz ] == npix ) {
7660 for ( inew = 1; inew <= ncur; inew++ ) {
7661 xcn = nmatx [ inew ];
7662 ycn = nmaty [ inew ];
7663 zcn = nmatz [ inew ];
7666 for ( jnew = 1; jnew <= 6; jnew++ ) {
7695 if ( y1 >= SYSIZE ) {
7697 }
else if ( y1 < 0 ) {
7702 if ( z1 >= SYSIZE ) {
7704 }
else if ( z1 < 0 ) {
7709 if ( ( x1 >= 0 ) && ( x1 < SYSIZE ) ) {
7711 px =
cx(x1, y1, z1, d1, d2, d3);
7712 py =
cy(x1, y1, z1, d1, d2, d3);
7713 pz =
cz(x1, y1, z1, d1, d2, d3);
7714 if (
mic [ px ] [ py ] [ pz ] == npix ) {
7719 printf(
"error in size of nnew \n");
7722 nnewx [ nnew ] = x1;
7723 nnewy [ nnew ] = y1;
7724 nnewz [ nnew ] = z1;
7733 for ( icur = 1; icur <= ncur; icur++ ) {
7734 nmatx [ icur ] = nnewx [ icur ];
7735 nmaty [ icur ] = nnewy [ icur ];
7736 nmatz [ icur ] = nnewz [ icur ];
7739 }
while ( nnew > 0 );
7745 for ( j1 = 0; j1 <
SYSIZE; j1++ ) {
7746 for ( k1 = 0; k1 <
SYSIZE; k1++ ) {
7747 px =
cx(xl, j1, k1, d1, d2, d3);
7748 py =
cy(xl, j1, k1, d1, d2, d3);
7749 pz =
cz(xl, j1, k1, d1, d2, d3);
7750 qx =
cx(xh, j1, k1, d1, d2, d3);
7751 qy =
cy(xh, j1, k1, d1, d2, d3);
7752 qz =
cz(xh, j1, k1, d1, d2, d3);
7753 if ( (
mic [ px ] [ py ] [ pz ] ==
BURNT ) && (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) ) {
7757 if (
mic [ px ] [ py ] [ pz ] ==
BURNT ) {
7758 mic [ px ] [ py ] [ pz ] =
BURNT + 1;
7761 if (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) {
7762 mic [ qx ] [ qy ] [ qz ] =
BURNT + 1;
7775 for ( i = 0; i <
SYSIZE; i++ ) {
7776 for ( j = 0; j <
SYSIZE; j++ ) {
7777 for ( k = 0; k <
SYSIZE; k++ ) {
7778 if (
mic [ i ] [ j ] [ k ] >=
BURNT ) {
7780 mic [ i ] [ j ] [ k ] = npix;
7781 }
else if (
mic [ i ] [ j ] [ k ] == npix ) {
7789 printf(
"Phase ID= %d \n", npix);
7790 printf(
"Number accessible from first surface = %ld \n", ntop);
7791 printf(
"Number contained in through pathways= %ld \n", nthrough);
7799 alpha_burn = 1. - ( mass_burn /
cemmass );
7806 if ( nthrough > 0 ) {
7818 ( void ) alpha_burn;
7831 long int ntop, nthrough, icur, inew, ncur, nnew, ntot, count_solid;
7832 int i, j, k, setyet;
7834 int *nmatx, *nmaty, *nmatz;
7835 int xl, xh, j1, k1, px, py, pz, qx, qy, qz;
7836 int xcn, ycn, zcn, x1, y1, z1, igood;
7838 int *nnewx, *nnewy, *nnewz;
7840 float mass_burn = 0.0, alpha_burn = 0.0, con_frac;
7857 for ( k = 0; k <
SYSIZE; k++ ) {
7858 for ( j = 0; j <
SYSIZE; j++ ) {
7859 for ( i = 0; i <
SYSIZE; i++ ) {
7860 newmat [ i ] [ j ] [ k ] =
mic [ i ] [ j ] [ k ];
7870 for ( k = 0; k <
SYSIZE; k++ ) {
7871 for ( j = 0; j <
SYSIZE; j++ ) {
7876 px =
cx(i, j, k, d1, d2, d3);
7877 py =
cy(i, j, k, d1, d2, d3);
7878 pz =
cz(i, j, k, d1, d2, d3);
7881 if ( (
mic [ px ] [ py ] [ pz ] ==
C3S ) ||
7882 (
mic [ px ] [ py ] [ pz ] ==
C2S ) ||
7883 (
mic [ px ] [ py ] [ pz ] ==
SLAG ) ||
7884 (
mic [ px ] [ py ] [ pz ] ==
ASG ) ||
7885 (
mic [ px ] [ py ] [ pz ] ==
CAS2 ) ||
7886 (
mic [ px ] [ py ] [ pz ] ==
POZZ ) ||
7887 (
mic [ px ] [ py ] [ pz ] ==
CSH ) ||
7888 (
mic [ px ] [ py ] [ pz ] ==
C3AH6 ) ||
7889 (
mic [ px ] [ py ] [ pz ] ==
ETTR ) ||
7891 (
mic [ px ] [ py ] [ pz ] ==
C3A ) ||
7892 (
mic [ px ] [ py ] [ pz ] ==
C4AF ) ) {
7905 for ( inew = 1; inew <= ncur; inew++ ) {
7906 xcn = nmatx [ inew ];
7907 ycn = nmaty [ inew ];
7908 zcn = nmatz [ inew ];
7910 qx =
cx(xcn, ycn, zcn, d1, d2, d3);
7911 qy =
cy(xcn, ycn, zcn, d1, d2, d3);
7912 qz =
cz(xcn, ycn, zcn, d1, d2, d3);
7915 for ( jnew = 1; jnew <= 6; jnew++ ) {
7944 if ( y1 >= SYSIZE ) {
7946 }
else if ( y1 < 0 ) {
7951 if ( z1 >= SYSIZE ) {
7953 }
else if ( z1 < 0 ) {
7958 if ( ( x1 >= 0 ) && ( x1 < SYSIZE ) ) {
7959 px =
cx(x1, y1, z1, d1, d2, d3);
7960 py =
cy(x1, y1, z1, d1, d2, d3);
7961 pz =
cz(x1, y1, z1, d1, d2, d3);
7964 if ( (
mic [ px ] [ py ] [ pz ] ==
CSH ) || (
mic [ px ] [ py ] [ pz ] ==
ETTRC4AF ) || (
mic [ px ] [ py ] [ pz ] ==
C3AH6 ) || (
mic [ px ] [ py ] [ pz ] ==
ETTR ) ) {
7969 printf(
"error in size of nnew %ld\n", nnew);
7972 nnewx [ nnew ] = x1;
7973 nnewy [ nnew ] = y1;
7974 nnewz [ nnew ] = z1;
7977 else if ( ( ( newmat [ qx ] [ qy ] [ qz ] ==
CSH ) || ( newmat [ qx ] [ qy ] [ qz ] ==
ETTRC4AF ) || ( newmat [ qx ] [ qy ] [ qz ] ==
C3AH6 ) || ( newmat [ qx ] [ qy ] [ qz ] ==
ETTR ) ) &&
7978 ( (
mic [ px ] [ py ] [ pz ] ==
C3S ) ||
7979 (
mic [ px ] [ py ] [ pz ] ==
C2S ) ||
7980 (
mic [ px ] [ py ] [ pz ] ==
CAS2 ) ||
7981 (
mic [ px ] [ py ] [ pz ] ==
SLAG ) ||
7982 (
mic [ px ] [ py ] [ pz ] ==
POZZ ) ||
7983 (
mic [ px ] [ py ] [ pz ] ==
ASG ) ||
7984 (
mic [ px ] [ py ] [ pz ] ==
C3A ) ||
7985 (
mic [ px ] [ py ] [ pz ] ==
C4AF ) ) ) {
7990 printf(
"error in size of nnew %ld\n", nnew);
7993 nnewx [ nnew ] = x1;
7994 nnewy [ nnew ] = y1;
7995 nnewz [ nnew ] = z1;
8000 else if ( (
micpart [ qx ] [ qy ] [ qz ] ==
micpart [ px ] [ py ] [ pz ] ) &&
8001 ( (
mic [ px ] [ py ] [ pz ] ==
C3S ) ||
8002 (
mic [ px ] [ py ] [ pz ] ==
C2S ) ||
8003 (
mic [ px ] [ py ] [ pz ] ==
POZZ ) ||
8004 (
mic [ px ] [ py ] [ pz ] ==
SLAG ) ||
8005 (
mic [ px ] [ py ] [ pz ] ==
ASG ) ||
8006 (
mic [ px ] [ py ] [ pz ] ==
CAS2 ) ||
8007 (
mic [ px ] [ py ] [ pz ] ==
C3A ) ||
8008 (
mic [ px ] [ py ] [ pz ] ==
C4AF ) ) && ( ( newmat [ qx ] [ qy ] [ qz ] ==
C3S ) ||
8009 ( newmat [ qx ] [ qy ] [ qz ] ==
C2S ) ||
8010 ( newmat [ qx ] [ qy ] [ qz ] ==
SLAG ) ||
8011 ( newmat [ qx ] [ qy ] [ qz ] ==
ASG ) ||
8012 ( newmat [ qx ] [ qy ] [ qz ] ==
POZZ ) ||
8013 ( newmat [ qx ] [ qy ] [ qz ] ==
CAS2 ) ||
8014 ( newmat [ qx ] [ qy ] [ qz ] ==
C3A ) ||
8015 ( newmat [ qx ] [ qy ] [ qz ] ==
C4AF ) ) ) {
8020 printf(
"error in size of nnew %ld\n", nnew);
8023 nnewx [ nnew ] = x1;
8024 nnewy [ nnew ] = y1;
8025 nnewz [ nnew ] = z1;
8039 for ( icur = 1; icur <= ncur; icur++ ) {
8040 nmatx [ icur ] = nnewx [ icur ];
8041 nmaty [ icur ] = nnewy [ icur ];
8042 nmatz [ icur ] = nnewz [ icur ];
8045 }
while ( nnew > 0 );
8051 for ( j1 = 0; j1 <
SYSIZE; j1++ ) {
8052 for ( k1 = 0; k1 <
SYSIZE; k1++ ) {
8053 px =
cx(xl, j1, k1, d1, d2, d3);
8054 py =
cy(xl, j1, k1, d1, d2, d3);
8055 pz =
cz(xl, j1, k1, d1, d2, d3);
8056 qx =
cx(xh, j1, k1, d1, d2, d3);
8057 qy =
cy(xh, j1, k1, d1, d2, d3);
8058 qz =
cz(xh, j1, k1, d1, d2, d3);
8059 if ( (
mic [ px ] [ py ] [ pz ] ==
BURNT ) && (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) ) {
8063 if (
mic [ px ] [ py ] [ pz ] ==
BURNT ) {
8064 mic [ px ] [ py ] [ pz ] =
BURNT + 1;
8067 if (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) {
8068 mic [ qx ] [ qy ] [ qz ] =
BURNT + 1;
8081 printf(
"Phase ID= Solid Phases \n");
8082 printf(
"Number accessible from first surface = %ld \n", ntop);
8083 printf(
"Number contained in through pathways= %ld \n", nthrough);
8090 alpha_burn = 1. - ( mass_burn /
cemmass );
8095 if ( count_solid > 0 ) {
8096 con_frac = ( float ) nthrough / (
float ) count_solid;
8104 if ( con_frac > 0.985 ) {
8111 for ( i = 0; i <
SYSIZE; i++ ) {
8112 for ( j = 0; j <
SYSIZE; j++ ) {
8113 for ( k = 0; k <
SYSIZE; k++ ) {
8114 if (
mic [ i ] [ j ] [ k ] >=
BURNT ) {
8115 mic [ i ] [ j ] [ k ] = newmat [ i ] [ j ] [ k ];
8133 ( void ) alpha_burn;
8141 int norig [ 100000 ], nleft [ 100000 ];
8143 char valmic, valmicorig;
8144 int valpart, partmax;
8149 for ( ix = 0; ix < 100000; ix++ ) {
8150 nleft [ ix ] = norig [ ix ] = 0;
8153 phydfile = fopen(
phrname,
"a");
8158 for ( ix = 0; ix <
SYSIZE; ix++ ) {
8159 for ( iy = 0; iy <
SYSIZE; iy++ ) {
8160 for ( iz = 0; iz <
SYSIZE; iz++ ) {
8161 if (
micpart [ ix ] [ iy ] [ iz ] != 0 ) {
8163 if ( valpart > partmax ) {
8167 valmic =
mic [ ix ] [
iy ] [ iz ];
8168 if ( ( valmic ==
C3S ) || ( valmic ==
C2S ) || ( valmic ==
C3A ) || ( valmic ==
C4AF ) ) {
8169 nleft [ valpart ] += 1;
8172 valmicorig =
micorig [ ix ] [
iy ] [ iz ];
8173 if ( ( valmicorig ==
C3S ) || ( valmicorig ==
C2S ) || ( valmicorig ==
C3A ) || ( valmicorig ==
C4AF ) ) {
8174 norig [ valpart ] += 1;
8182 for ( ix = 100; ix <= partmax; ix++ ) {
8184 if ( norig [ ix ] != 0 ) {
8185 alpart = 1. - ( float ) nleft [ ix ] / (
float ) norig [ ix ];
8188 fprintf(phydfile,
"%d %d %d %.3f\n", ix, norig [ ix ], nleft [ ix ], alpart);
8199 int plok, sumnew, xl1, yl1, zl1, act1;
8210 plok = ( int ) ( 6. *
ran1(
seed) );
8211 if ( ( plok > 5 ) || ( plok < 0 ) ) {
8223 if ( sumold % 2 != 0 ) {
8235 if ( sumold % 3 != 0 ) {
8247 if ( sumold % 5 != 0 ) {
8259 if ( sumold % 7 != 0 ) {
8271 if ( sumold % 11 != 0 ) {
8283 if ( sumold % 13 != 0 ) {
8310 int ixe, iye, ize, edgeback, x2, y2, z2, check;
8317 for ( ixe = ( -1 ); ixe <= 1; ixe++ ) {
8319 for ( iye = ( -1 ); iye <= 1; iye++ ) {
8321 for ( ize = ( -1 ); ize <= 1; ize++ ) {
8322 if ( ( ixe != 0 ) || ( iye != 0 ) || ( ize != 0 ) ) {
8327 }
else if ( x2 >=
SYSIZE ) {
8333 }
else if ( y2 >=
SYSIZE ) {
8339 }
else if ( z2 >=
SYSIZE ) {
8343 check =
mic [ x2 ] [ y2 ] [ z2 ];
8344 if ( ( check != ph1 ) && ( check != ph2 ) && ( check != ph3 ) ) {
8353 return ( edgeback );
8361 int numnear, xchr, ychr, zchr, fchr, check, msface;
8368 while ( fchr == 0 ) {
8386 check =
mic [ xchr ] [ ychr ] [ zchr ];
8393 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8394 mic [ xchr ] [ ychr ] [ zchr ] =
CSH;
8399 msface = ( int ) ( 3. *
ran1(
seed) + 1. );
8404 faces [ xchr ] [ ychr ] [ zchr ] = msface;
8422 int xnew, ynew, znew, action, sumback, sumin, check;
8423 int msface, mstest = 0, mstest2 = 0;
8424 float prcsh, prcsh1, prcsh2, prtest;
8432 sumback =
moveone(& xnew, & ynew, & znew, & action, sumin);
8435 if ( xnew != xcur ) {
8440 if ( ynew != ycur ) {
8445 if ( znew != zcur ) {
8451 if ( action == 0 ) {
8452 printf(
"Error in value of action \n");
8455 check =
mic [ xnew ] [ ynew ] [ znew ];
8461 if ( ( check ==
CSH ) && ( (
cshgeom == 0 ) || (
faces [ xnew ] [ ynew ] [ znew ] == 0 ) || (
faces [ xnew ] [ ynew ] [ znew ] == mstest ) || (
faces [ xnew ] [ ynew ] [ znew ] == mstest2 ) ) ) {
8467 if ( prcsh1 <= prtest ) {
8468 mic [ xcur ] [ ycur ] [ zcur ] =
CSH;
8470 faces [ xcur ] [ ycur ] [ zcur ] =
faces [ xnew ] [ ynew ] [ znew ];
8482 if ( prtest > 1.0 ) {
8484 if ( prcsh2 < ( prtest - 1.0 ) ) {
8492 else if ( ( check ==
SLAGCSH ) || ( check ==
POZZCSH ) || ( finalstep == 1 ) ||
8493 ( ( ( check ==
C3S ) || ( check ==
C2S ) ) && ( prcsh < 0.001 ) ) ||
8494 ( ( ( check ==
C3A ) || ( check ==
C4AF ) ) && ( prcsh < 0.2 ) ) ||
8495 ( ( check ==
CH ) && ( prcsh < 0.01 ) ) ||
8502 if ( prcsh1 <= prtest ) {
8503 mic [ xcur ] [ ycur ] [ zcur ] =
CSH;
8506 msface = ( int ) ( 2. *
ran1(
seed) + 1. );
8511 if ( msface == 1 ) {
8512 faces [ xcur ] [ ycur ] [ zcur ] = mstest;
8514 faces [ xcur ] [ ycur ] [ zcur ] = mstest2;
8527 if ( prtest > 1.0 ) {
8529 if ( prcsh2 < ( prtest - 1.0 ) ) {
8537 if ( action != 0 ) {
8560 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
8569 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
8575 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
8576 if ( newact == 0 ) {
8577 printf(
"Error in value of newact in extfh3 \n");
8580 check =
mic [ xchr ] [ ychr ] [ zchr ];
8585 mic [ xchr ] [ ychr ] [ zchr ] =
FH3;
8597 while ( fchr == 0 ) {
8615 check =
mic [ xchr ] [ ychr ] [ zchr ];
8621 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8622 mic [ xchr ] [ ychr ] [ zchr ] =
FH3;
8641 int check, newact, multf, numnear, sump, xchr, ychr, zchr, fchr, i1;
8642 int numalum, numsil;
8643 float pneigh, ptest;
8653 for ( i1 = 1; ( ( i1 <= 1000 ) && ( fchr == 0 ) ); i1++ ) {
8659 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
8660 if ( newact == 0 ) {
8661 printf(
"Error in value of action \n");
8664 check =
mic [ xchr ] [ ychr ] [ zchr ];
8672 numsil = 26 - numsil;
8676 numalum = 26 - numalum;
8680 numalum = 26 - numalum;
8683 pneigh = ( float ) ( numnear + 1 ) / 26.0;
8685 if ( numalum <= 1 ) {
8689 if ( numalum >= 2 ) {
8693 if ( numalum >= 3 ) {
8697 if ( numalum >= 5 ) {
8703 if ( pneigh >= ptest ) {
8705 mic [ xchr ] [ ychr ] [ zchr ] =
ETTR;
8723 while ( fchr == 0 ) {
8742 check =
mic [ xchr ] [ ychr ] [ zchr ];
8746 numsil = 26 - numsil;
8755 if ( ( tries > 5000 ) || ( ( numnear < 26 ) && ( numsil < 1 ) ) ) {
8757 mic [ xchr ] [ ychr ] [ zchr ] =
ETTR;
8781 int numnear, xchr, ychr, zchr, fchr, check;
8788 while ( fchr == 0 ) {
8806 check =
mic [ xchr ] [ ychr ] [ zchr ];
8813 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8814 mic [ xchr ] [ ychr ] [ zchr ] =
CH;
8828 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
8837 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
8843 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
8844 if ( newact == 0 ) {
8845 printf(
"Error in value of newact in extfh3 \n");
8848 check =
mic [ xchr ] [ ychr ] [ zchr ];
8865 while ( fchr == 0 ) {
8883 check =
mic [ xchr ] [ ychr ] [ zchr ];
8889 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8907 int xnew, ynew, znew, action, sumback, sumin, check = -1;
8908 int nexp, iexp, xexp, yexp, zexp, newact, ettrtype;
8909 float pgen, pexp, pext, p2diff;
8915 if ( ( nucprgyp >= pgen ) || ( finalstep == 1 ) ) {
8930 sumback =
moveone(& xnew, & ynew, & znew, & action, sumin);
8932 if ( action == 0 ) {
8933 printf(
"Error in value of action \n");
8936 check =
mic [ xnew ] [ ynew ] [ znew ];
8958 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
8966 count [ check ] -= 1;
8973 if ( pexp <= 0.569 ) {
8974 if ( ettrtype == 0 ) {
8975 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
8986 if ( check ==
C3A ) {
8987 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
8990 if ( ettrtype == 0 ) {
9008 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9009 newact =
extettr(xexp, yexp, zexp, ettrtype);
9061 if ( pexp <= 0.6935 ) {
9062 newact =
extettr(xexp, yexp, zexp, ettrtype);
9078 if ( pexp <= 0.8174 ) {
9085 if ( pext < 0.2584 ) {
9091 if ( pext < 0.5453 ) {
9092 extfh3(xnew, ynew, znew);
9097 mic [ xnew ] [ ynew ] [ znew ] =
C4AF;
9107 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9108 newact =
extettr(xexp, yexp, zexp, 1);
9160 if ( pexp <= 0.6935 ) {
9161 newact =
extettr(xexp, yexp, zexp, 1);
9168 if ( action != 0 ) {
9193 int xnew, ynew, znew, action, sumback, sumin, check = -1;
9194 int nexp, iexp, xexp, yexp, zexp, newact, ettrtype;
9195 float pgen, pexp, pext, p2diff;
9201 if ( ( nucprgyp >= pgen ) || ( finalstep == 1 ) ) {
9217 sumback =
moveone(& xnew, & ynew, & znew, & action, sumin);
9219 if ( action == 0 ) {
9220 printf(
"Error in value of action \n");
9223 check =
mic [ xnew ] [ ynew ] [ znew ];
9245 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
9253 count [ check ] -= 1;
9260 if ( pexp <= 0.5583 ) {
9261 if ( ettrtype == 0 ) {
9262 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
9273 if ( check ==
C3A ) {
9274 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
9277 if ( ettrtype == 0 ) {
9295 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9296 newact =
extettr(xexp, yexp, zexp, ettrtype);
9348 if ( pexp <= 0.6053 ) {
9349 newact =
extettr(xexp, yexp, zexp, ettrtype);
9365 if ( pexp <= 0.802 ) {
9372 if ( pext < 0.2584 ) {
9378 if ( pext < 0.5453 ) {
9379 extfh3(xnew, ynew, znew);
9384 mic [ xnew ] [ ynew ] [ znew ] =
C4AF;
9394 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9395 newact =
extettr(xexp, yexp, zexp, 1);
9447 if ( pexp <= 0.6053 ) {
9448 newact =
extettr(xexp, yexp, zexp, 1);
9455 if ( action != 0 ) {
9478 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
9487 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
9493 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
9494 if ( newact == 0 ) {
9495 printf(
"Error in value of newact in extfreidel \n");
9498 check =
mic [ xchr ] [ ychr ] [ zchr ];
9515 while ( fchr == 0 ) {
9534 check =
mic [ xchr ] [ ychr ] [ zchr ];
9540 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
9559 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
9568 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
9574 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
9575 if ( newact == 0 ) {
9576 printf(
"Error in value of newact in extstrat \n");
9579 check =
mic [ xchr ] [ ychr ] [ zchr ];
9584 mic [ xchr ] [ ychr ] [ zchr ] =
STRAT;
9596 while ( fchr == 0 ) {
9615 check =
mic [ xchr ] [ ychr ] [ zchr ];
9621 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
9622 mic [ xchr ] [ ychr ] [ zchr ] =
STRAT;
9640 int check, xnew, ynew, znew, action, nexp, iexp;
9641 int xexp, yexp, zexp, newact, sumold, sumgarb, ettrtype;
9642 float pexp, pext, p2diff;
9648 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFGYP ) {
9658 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
9659 if ( action == 0 ) {
9660 printf(
"Error in value of action in movegyp \n");
9663 check =
mic [ xnew ] [ ynew ] [ znew ];
9672 mic [ xcur ] [ ycur ] [ zcur ] =
ABSGYP;
9683 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
9691 count [ check ] -= 1;
9698 if ( pexp <= 0.40 ) {
9699 if ( ettrtype == 0 ) {
9700 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
9711 if ( check ==
C3A ) {
9712 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
9715 if ( ettrtype == 0 ) {
9733 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9734 newact =
extettr(xexp, yexp, zexp, ettrtype);
9786 if ( pexp <= 0.30 ) {
9787 newact =
extettr(xexp, yexp, zexp, ettrtype);
9803 if ( pexp <= 0.575 ) {
9810 if ( pext < 0.2584 ) {
9816 if ( pext < 0.5453 ) {
9817 extfh3(xnew, ynew, znew);
9822 mic [ xnew ] [ ynew ] [ znew ] =
C4AF;
9832 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9833 newact =
extettr(xexp, yexp, zexp, 1);
9885 if ( pexp <= 0.30 ) {
9886 newact =
extettr(xexp, yexp, zexp, 1);
9894 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
9898 mic [ xcur ] [ ycur ] [ zcur ] =
GYPSUM;
9901 if ( action != 0 ) {
9925 int check, xnew, ynew, znew, action, nexp, iexp;
9926 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
9944 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
9945 if ( action == 0 ) {
9946 printf(
"Error in value of action in movecacl2 \n");
9949 check =
mic [ xnew ] [ ynew ] [ znew ];
9958 count [ check ] -= 1;
9965 if ( pexp <= 0.5793 ) {
9980 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10033 if ( pexp <= 0.3295 ) {
10039 else if ( check ==
C4AF ) {
10051 if ( pexp <= 0.4033 ) {
10058 if ( pext < 0.6412 ) {
10064 if ( pext < 0.3522 ) {
10065 extfh3(xnew, ynew, znew);
10068 extfh3(xnew, ynew, znew);
10079 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10082 switch ( newact ) {
10132 if ( pexp <= 0.3176 ) {
10141 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10145 mic [ xcur ] [ ycur ] [ zcur ] =
CACL2;
10148 if ( action != 0 ) {
10176 int check, xnew, ynew, znew, action, nexp, iexp;
10177 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
10185 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFCAS2 ) {
10195 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10196 if ( action == 0 ) {
10197 printf(
"Error in value of action in movecas2 \n");
10200 check =
mic [ xnew ] [ ynew ] [ znew ];
10207 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
10216 if ( pexp <= 0.886 ) {
10217 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
10219 count [ check ] -= 1;
10228 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10229 newact =
extstrat(xexp, yexp, zexp);
10231 switch ( newact ) {
10281 if ( pexp <= 0.286 ) {
10282 newact =
extstrat(xexp, yexp, zexp);
10287 else if ( check ==
C4AF ) {
10288 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
10299 if ( pexp <= 0.786 ) {
10300 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
10307 if ( pext < 0.329 ) {
10314 if ( pext < 0.6938 ) {
10315 extfh3(xnew, ynew, znew);
10327 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10328 newact =
extstrat(xexp, yexp, zexp);
10330 switch ( newact ) {
10380 if ( pexp <= 0.37 ) {
10381 newact =
extstrat(xexp, yexp, zexp);
10389 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10393 mic [ xcur ] [ ycur ] [ zcur ] =
CAS2;
10396 if ( action != 0 ) {
10424 int check, xnew, ynew, znew, action, nexp, iexp;
10425 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
10433 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFAS ) {
10443 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10444 if ( action == 0 ) {
10445 printf(
"Error in value of action in moveas \n");
10448 check =
mic [ xnew ] [ ynew ] [ znew ];
10452 if ( ( check ==
CH ) || ( check ==
DIFFCH ) ) {
10455 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
10457 count [ check ] -= 1;
10464 if ( pexp <= 0.7538 ) {
10465 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
10479 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10480 newact =
extstrat(xexp, yexp, zexp);
10482 switch ( newact ) {
10532 if ( pexp <= 0.326 ) {
10533 newact =
extstrat(xexp, yexp, zexp);
10539 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10543 mic [ xcur ] [ ycur ] [ zcur ] =
ASG;
10546 if ( action != 0 ) {
10550 mic [ xnew ] [ ynew ] [ znew ] =
DIFFAS;
10574 int check, xnew, ynew, znew, action;
10575 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
10593 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10594 if ( action == 0 ) {
10595 printf(
"Error in value of action in moveas \n");
10598 check =
mic [ xnew ] [ ynew ] [ znew ];
10603 if ( check ==
AFM ) {
10607 if ( pexp <= 0.479192 ) {
10608 mic [ xnew ] [ ynew ] [ znew ] =
AFMC;
10611 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
10615 count [ check ] -= 1;
10621 if ( pexp <= 0.078658 ) {
10622 mic [ xcur ] [ ycur ] [ zcur ] =
AFMC;
10637 if ( pexp <= 0.26194 ) {
10638 newact =
extettr(xexp, yexp, zexp, 0);
10644 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10648 mic [ xcur ] [ ycur ] [ zcur ] =
CACO3;
10651 if ( action != 0 ) {
10679 int check, sump, xchr, ychr, zchr, fchr, i1, newact, numnear;
10688 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
10694 sump *=
moveone(& xchr, & ychr, & zchr, & newact, sump);
10695 if ( newact == 0 ) {
10696 printf(
"Error in value of newact in extafm \n");
10699 check =
mic [ xchr ] [ ychr ] [ zchr ];
10703 mic [ xchr ] [ ychr ] [ zchr ] =
AFM;
10713 while ( fchr == 0 ) {
10731 check =
mic [ xchr ] [ ychr ] [ zchr ];
10738 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
10739 mic [ xchr ] [ ychr ] [ zchr ] =
AFM;
10754 int check, xnew, ynew, znew, action;
10755 int sumold, sumgarb;
10756 float pexp, pafm, pgrow;
10760 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFETTR ) {
10771 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10772 if ( action == 0 ) {
10773 printf(
"Error in value of action in moveettr \n");
10776 check =
mic [ xnew ] [ ynew ] [ znew ];
10780 if ( check ==
C4AF ) {
10782 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
10792 if ( pexp <= 0.278 ) {
10793 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
10798 if ( pafm <= 0.3241 ) {
10804 if ( pafm <= 0.4313 ) {
10805 extfh3(xnew, ynew, znew);
10807 }
else if ( pexp <= 0.348 ) {
10808 mic [ xnew ] [ ynew ] [ znew ] =
FH3;
10817 else if ( ( check ==
C3A ) || ( check ==
DIFFC3A ) ) {
10820 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
10823 count [ check ] -= 1;
10829 if ( pexp <= 0.2424 ) {
10830 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
10836 if ( check ==
C3A ) {
10837 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
10850 if ( pexp <= pafm ) {
10851 extafm(xcur, ycur, zcur);
10855 else if ( check ==
ETTR ) {
10858 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
10867 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10871 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
10874 if ( action != 0 ) {
10897 int check, sump, xchr, ychr, zchr, fchr, i1, action, numnear;
10906 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
10912 sump *=
moveone(& xchr, & ychr, & zchr, & action, sump);
10913 if ( action == 0 ) {
10914 printf(
"Error in value of action in extpozz \n");
10917 check =
mic [ xchr ] [ ychr ] [ zchr ];
10931 while ( fchr == 0 ) {
10949 check =
mic [ xchr ] [ ychr ] [ zchr ];
10955 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
10971 int check, xnew, ynew, znew, action, sumold, sumgarb;
10977 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
10979 mic [ xcur ] [ ycur ] [ zcur ] =
FH3;
10989 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10990 if ( action == 0 ) {
10991 printf(
"Error in value of action in movefh3 \n");
10994 check =
mic [ xnew ] [ ynew ] [ znew ];
10997 if ( check ==
FH3 ) {
10998 mic [ xcur ] [ ycur ] [ zcur ] =
FH3;
11004 if ( action != 0 ) {
11028 int check, xnew, ynew, znew, action, sumgarb, sumold;
11029 float pexp, pgen, pfix;
11033 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
11035 mic [ xcur ] [ ycur ] [ zcur ] =
CH;
11045 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
11046 if ( action == 0 ) {
11047 printf(
"Error in value of action in movech \n");
11050 check =
mic [ xnew ] [ ynew ] [ znew ];
11053 if ( ( check ==
CH ) && ( pgen <=
CHGROW ) ) {
11054 mic [ xcur ] [ ycur ] [ zcur ] =
CH;
11062 mic [ xcur ] [ ycur ] [ zcur ] =
CH;
11069 else if ( ( pgen <=
ppozz ) && ( check ==
POZZ ) && (
npr <= (
int ) ( (
float )
nfill * 1.35 ) ) ) {
11079 if ( pfix <= ( 1. / 1.35 ) ) {
11091 if ( pexp <= 0.05466 ) {
11094 }
else if ( check ==
DIFFAS ) {
11096 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
11104 if ( pfix <= 0.7538 ) {
11105 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
11114 if ( pexp <= 0.5035 ) {
11119 if ( action != 0 ) {
11123 mic [ xnew ] [ ynew ] [ znew ] =
DIFFCH;
11143 int check, sump, xchr, ychr, zchr, fchr, i1, action, numnear;
11152 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
11158 sump *=
moveone(& xchr, & ychr, & zchr, & action, sump);
11159 if ( action == 0 ) {
11160 printf(
"Error in action value in extc3ah6 \n");
11163 check =
mic [ xchr ] [ ychr ] [ zchr ];
11167 mic [ xchr ] [ ychr ] [ zchr ] =
C3AH6;
11176 while ( fchr == 0 ) {
11193 check =
mic [ xchr ] [ ychr ] [ zchr ];
11199 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
11200 mic [ xchr ] [ ychr ] [ zchr ] =
C3AH6;
11215 int check, xnew, ynew, znew, action, sumgarb, sumold;
11216 int xexp, yexp, zexp, nexp, iexp, newact;
11217 float pgen, pexp, pafm, pgrow, p2diff;
11220 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFC3A ) {
11229 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
11231 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11238 if ( pexp <= 0.69 ) {
11248 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
11249 if ( action == 0 ) {
11250 printf(
"Error in value of action in movec3a \n");
11253 check =
mic [ xnew ] [ ynew ] [ znew ];
11256 if ( check ==
C3AH6 ) {
11261 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11268 if ( pexp <= 0.69 ) {
11277 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
11286 if ( pexp <= 0.40 ) {
11287 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
11303 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11304 newact =
extettr(xexp, yexp, zexp, 0);
11306 switch ( newact ) {
11356 if ( pexp <= 0.30 ) {
11357 newact =
extettr(xexp, yexp, zexp, 0);
11364 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
11373 if ( pexp <= 0.5583 ) {
11374 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
11390 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11391 newact =
extettr(xexp, yexp, zexp, 0);
11393 switch ( newact ) {
11443 if ( pexp <= 0.6053 ) {
11444 newact =
extettr(xexp, yexp, zexp, 0);
11451 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
11460 if ( pexp <= 0.569 ) {
11461 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
11477 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11478 newact =
extettr(xexp, yexp, zexp, 0);
11480 switch ( newact ) {
11530 if ( pexp <= 0.6935 ) {
11531 newact =
extettr(xexp, yexp, zexp, 0);
11547 if ( pexp <= 0.5793 ) {
11562 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11565 switch ( newact ) {
11615 if ( pexp <= 0.3295 ) {
11623 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
11632 if ( pexp <= 0.886 ) {
11633 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
11648 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11649 newact =
extstrat(xexp, yexp, zexp);
11651 switch ( newact ) {
11701 if ( pexp <= 0.286 ) {
11702 newact =
extstrat(xexp, yexp, zexp);
11713 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
11716 count [ check ] -= 1;
11721 if ( pexp <= 0.2424 ) {
11722 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
11733 if ( pexp <= pafm ) {
11734 extafm(xnew, ynew, znew);
11738 if ( ( action != 0 ) && ( action != 7 ) ) {
11762 int check, xnew, ynew, znew, action, sumgarb, sumold;
11763 int xexp, yexp, zexp, nexp, iexp, newact;
11764 float pgen, pexp, pafm, pgrow, p2diff;
11767 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFC4A ) {
11776 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
11778 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11785 if ( pexp <= 0.69 ) {
11795 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
11796 if ( action == 0 ) {
11797 printf(
"Error in value of action in movec4a \n");
11800 check =
mic [ xnew ] [ ynew ] [ znew ];
11803 if ( check ==
C3AH6 ) {
11808 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11815 if ( pexp <= 0.69 ) {
11833 if ( pexp <= 0.40 ) {
11850 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11851 newact =
extettr(xexp, yexp, zexp, 1);
11853 switch ( newact ) {
11903 if ( pexp <= 0.30 ) {
11904 newact =
extettr(xexp, yexp, zexp, 1);
11920 if ( pexp <= 0.5583 ) {
11937 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11938 newact =
extettr(xexp, yexp, zexp, 1);
11940 switch ( newact ) {
11990 if ( pexp <= 0.6053 ) {
11991 newact =
extettr(xexp, yexp, zexp, 1);
12007 if ( pexp <= 0.569 ) {
12024 for ( iexp = 1; iexp <= nexp; iexp++ ) {
12025 newact =
extettr(xexp, yexp, zexp, 1);
12027 switch ( newact ) {
12077 if ( pexp <= 0.6935 ) {
12078 newact =
extettr(xexp, yexp, zexp, 1);
12094 if ( pexp <= 0.5793 ) {
12109 for ( iexp = 1; iexp <= nexp; iexp++ ) {
12112 switch ( newact ) {
12162 if ( pexp <= 0.3295 ) {
12170 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
12179 if ( pexp <= 0.886 ) {
12180 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
12195 for ( iexp = 1; iexp <= nexp; iexp++ ) {
12196 newact =
extstrat(xexp, yexp, zexp);
12198 switch ( newact ) {
12248 if ( pexp <= 0.286 ) {
12249 newact =
extstrat(xexp, yexp, zexp);
12260 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
12263 count [ check ] -= 1;
12268 if ( pexp <= 0.2424 ) {
12269 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
12280 if ( pexp <= pafm ) {
12281 extafm(xnew, ynew, znew);
12285 if ( ( action != 0 ) && ( action != 7 ) ) {
12307 void CemhydMatStatus :: hydrate(
int fincyc,
int stepmax,
float chpar1,
float chpar2,
float hgpar1,
float hgpar2,
float fhpar1,
float fhpar2,
float gypar1,
float gypar2)
12309 int xpl, ypl, zpl, phpl, agepl, xpnew, ypnew, zpnew;
12310 float chprob, c3ah6prob, fh3prob, gypprob;
12311 long int nleft, ntodo, ndale;
12312 int istep, termflag, reactf = -1;
12314 struct ants *curant, *antgone;
12321 for ( istep = 1; ( ( istep <= stepmax ) && ( nleft > 0 ) ); istep++ ) {
12322 if ( ( fincyc == 1 ) && ( istep == stepmax ) ) {
12331 chprob = chpar1 * ( 1. - beterm );
12333 c3ah6prob = hgpar1 * ( 1. - beterm );
12335 fh3prob = fhpar1 * ( 1. - beterm );
12337 gypprob = gypar1 * ( 1. - beterm );
12341 while ( curant != NULL ) {
12347 agepl = curant->cycbirth;
12353 reactf =
movecsh(xpl, ypl, zpl, termflag, agepl);
12354 }
else if ( phpl ==
DIFFANH ) {
12357 reactf =
moveanh(xpl, ypl, zpl, termflag, gypprob);
12358 }
else if ( phpl == DIFFHEM ) {
12361 reactf =
movehem(xpl, ypl, zpl, termflag, gypprob);
12362 }
else if ( phpl == DIFFCH ) {
12365 reactf =
movech(xpl, ypl, zpl, termflag, chprob);
12366 }
else if ( phpl == DIFFFH3 ) {
12369 reactf =
movefh3(xpl, ypl, zpl, termflag, fh3prob);
12370 }
else if ( phpl ==
DIFFGYP ) {
12373 reactf =
movegyp(xpl, ypl, zpl, termflag);
12374 }
else if ( phpl == DIFFC3A ) {
12377 reactf =
movec3a(xpl, ypl, zpl, termflag, c3ah6prob);
12378 }
else if ( phpl ==
DIFFC4A ) {
12381 reactf =
movec4a(xpl, ypl, zpl, termflag, c3ah6prob);
12385 reactf =
moveettr(xpl, ypl, zpl, termflag);
12389 reactf =
movecacl2(xpl, ypl, zpl, termflag);
12393 reactf =
movecas2(xpl, ypl, zpl, termflag);
12394 }
else if ( phpl ==
DIFFAS ) {
12397 reactf =
moveas(xpl, ypl, zpl, termflag);
12401 reactf =
movecaco3(xpl, ypl, zpl, termflag);
12403 printf(
"Error in ID of phase \n");
12407 if ( reactf != 0 ) {
12414 switch ( reactf ) {
12424 if ( xpnew >=
SYSIZE ) {
12438 if ( ypnew >=
SYSIZE ) {
12452 if ( zpnew >=
SYSIZE ) {
12466 curant = curant->nextant;
12470 if ( ndale == 1 ) {
12473 ( curant->prevant )->nextant = curant->
nextant;
12476 if ( curant->nextant != NULL ) {
12477 ( curant->nextant )->prevant = curant->
prevant;
12502 float err, dxold, cdx, abx;
12503 fcomplex_cem sq, h,
gp, gm, g2, g, b, d, dx, f, x1;
12506 for ( iter = 1; iter <=
MAXIT; iter++ ) {
12511 for ( j = m - 1; j >= 0; j-- ) {
12515 err =
Cabs(b) + abx * err;
12519 if (
Cabs(b) <= err ) {
12534 x1 =
Csub(* x, dx);
12535 if ( x->
r == x1.
r && x->
i == x1.
i ) {
12541 if ( iter > 6 && cdx >= dxold ) {
12547 if ( cdx <= eps *
Cabs(* x) ) {
12553 nrerror(
"Too many iterations in routine CemhydMat - LAGUER\n");
12566 for ( j = 0; j <= m; j++ ) {
12567 ad [ j ] = a [ j ];
12570 for ( j = m; j >= 1; j-- ) {
12573 if ( fabs(x.
i) <= ( 2.0 *
EPSP * fabs(x.
r) ) ) {
12579 for ( jj = j - 1; jj >= 0; jj-- ) {
12587 for ( j = 1; j <= m; j++ ) {
12592 for ( j = 2; j <= m; j++ ) {
12594 for ( i = j - 1; i >= 1; i-- ) {
12595 if ( roots [ i ].r <= x.
r ) {
12599 roots [ i + 1 ] = roots [ i ];
12602 roots [ i + 1 ] = x;
12612 int j, syngen_change = 0, syn_old = 0, counter = 0;
12613 double concnaplus, conckplus;
12614 double concohminus = 0., A, B, C, conctest, concsulfate1;
12615 double volpore, grams_cement;
12616 double releasedna, releasedk, activitySO4, activityK, test_precip;
12617 double activityCa, activityOH, Istrength, Anow, Bnow, Inew;
12618 double conductivity = 0.0;
12620 float sumbest, sumtest, pozzreact, KspCH;
12640 volpore /= grams_cement;
12642 pozzreact = ( ( float )
npr / 1.35 ) * MASSFACTOR * MASSFACTOR * MASSFACTOR *
specgrav [
POZZ ];
12647 releasedk /=
MMK2O;
12655 releasedk /=
MMK2O;
12668 activityCa = activityOH = activitySO4 = activityK = 1.0;
12670 if ( ( ( concnaplus + conckplus ) > 0.0 ) && (
soluble [
ETTR ] == 0 ) ) {
12673 if ( Istrength < 1.0 ) {
12677 while ( ( fabs(Istrength - Inew) / Istrength ) > 0.10 ) {
12679 if ( Istrength < 1.0 ) {
12686 activityCa = ( -Anow *
zCa *
zCa * sqrt(Istrength) ) / ( 1. +
aCa * Bnow * sqrt(Istrength) );
12687 activityCa += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zCa *
zCa * Istrength / sqrt(1000.);
12688 activityCa = exp(activityCa);
12689 activityOH = ( -Anow *
zOH *
zOH * sqrt(Istrength) ) / ( 1. +
aOH * Bnow * sqrt(Istrength) );
12690 activityOH += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zOH *
zOH * Istrength / sqrt(1000.);
12691 activityOH = exp(activityOH);
12692 activityK = ( -Anow *
zK *
zK * sqrt(Istrength) ) / ( 1. +
aK * Bnow * sqrt(Istrength) );
12693 activityK += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zK *
zK * Istrength / sqrt(1000.);
12694 activityK = exp(activityK);
12695 activitySO4 = ( -Anow *
zSO4 *
zSO4 * sqrt(Istrength) ) / ( 1. +
aSO4 * Bnow * sqrt(Istrength) );
12696 activitySO4 += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zSO4 *
zSO4 * Istrength / sqrt(1000.);
12697 activitySO4 = exp(activitySO4);
12702 A = ( -KspCH / ( activityCa * activityOH * activityOH ) );
12703 B = conckplus + concnaplus;
12704 C = ( -2. *
KspGypsum / ( activityCa * activitySO4 ) );
12705 concohminus = conckplus + concnaplus;
12720 zroots(coef, 4, roots, 1);
12723 for ( j = 1; j <= 4; j++ ) {
12726 if ( ( ( roots [ j ].i ) == 0.0 ) && ( ( roots [ j ].r ) > 0.0 ) ) {
12727 conctest = sqrt( KspCH / ( roots [ j ].r * activityCa * activityOH * activityOH ) );
12728 concsulfate1 =
KspGypsum / ( roots [ j ].
r * activityCa * activitySO4 );
12729 sumtest = concnaplus + conckplus + 2. * roots [ j ].
r - conctest - 2. * concsulfate1;
12730 if ( fabs(sumtest) < sumbest ) {
12731 sumbest = fabs(sumtest);
12732 concohminus = conctest;
12747 if ( Istrength < 1.0 ) {
12751 while ( ( fabs(Istrength - Inew) / Istrength ) > 0.10 && counter < 4 ) {
12757 activityCa = ( -Anow *
zCa *
zCa * sqrt(Istrength) ) / ( 1. +
aCa * Bnow * sqrt(Istrength) );
12758 activityCa += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zCa *
zCa * Istrength / sqrt(1000.);
12759 activityCa = exp(activityCa);
12760 activityOH = ( -Anow *
zOH *
zOH * sqrt(Istrength) ) / ( 1. +
aOH * Bnow * sqrt(Istrength) );
12761 activityOH += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zOH *
zOH * Istrength / sqrt(1000.);
12762 activityOH = exp(activityOH);
12763 activityK = ( -Anow *
zK *
zK * sqrt(Istrength) ) / ( 1. +
aK * Bnow * sqrt(Istrength) );
12764 activityK += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zK *
zK * Istrength / sqrt(1000.);
12765 activityK = exp(activityK);
12767 concohminus = conckplus + concnaplus;
12768 if ( (
conccaplus ) > ( 0.1 * ( concohminus ) ) ) {
12772 conccaplus = ( KspCH / ( activityCa * activityOH * activityOH * concohminus * concohminus ) );
12784 if ( syn_old != 2 ) {
12785 test_precip = conckplus * conckplus * activityK * activityK;
12790 printf(
"Syngenite precipitating at cycle %d\n",
icyc);
12792 syngen_change = syn_old = 1;
12794 if ( conckplus > 0.002 ) {
12795 conckplus -= 0.001;
12797 }
else if ( conckplus > 0.0002 ) {
12798 conckplus -= 0.0001;
12812 syngen_change = syn_old = 2;
12815 conckplus += 0.001 *
KperSyn;
12823 }
while ( syngen_change != 0 );
12825 if ( concohminus < ( 0.0000001 ) ) {
12826 concohminus = 0.0000001;
12827 conccaplus = ( KspCH / ( activityCa * activityOH * activityOH * concohminus * concohminus ) );
12830 pH_cur = 14.0 + log10(concohminus * activityOH);
12833 Istrength /= 1000.;
12835 conductivity +=
zOH * concohminus * (
lambdaOH_0 / ( 1. +
GOH * sqrt(Istrength) ) );
12836 conductivity +=
zNa * concnaplus * (
lambdaNa_0 / ( 1. +
GNa * sqrt(Istrength) ) );
12837 conductivity +=
zK * conckplus * (
lambdaK_0 / ( 1. +
GK * sqrt(Istrength) ) );
12843 fprintf(
pHfile,
"%d %.4f %f %.4f %f %f %f %f %f %.4f %.4f %.4f %.4f %f\n",
cyccnt - 1,
time_cur,
alpha_cur,
pH_cur, conductivity, concnaplus, conckplus, conccaplus, concsulfate, activityCa, activityOH, activitySO4, activityK,
moles_syn_precip);
12851 if ( ( phase >=
C3S && phase <=
ABSGYP ) || ( phase ==
HDCSH ) ) {
12860 long int icur, inew, ncur, nnew;
12861 int i, j, k, cnt_perc, cnt_tot, kh;
12862 int *nmatx, *nmaty, *nmatz;
12863 int xl, xh, j1, k1, px, py, pz, qx, qy, qz;
12864 int xcn, ycn, zcn, x1, y1, z1, igood;
12865 int *nnewx, *nnewy, *nnewz;
12867 int phase_temp [ 51 ];
12888 for ( k = 0; k <
SYSIZE; k++ ) {
12889 for ( j = 0; j <
SYSIZE; j++ ) {
12890 for ( i = 0; i <
SYSIZE; i++ ) {
12891 newmat [ i ] [ j ] [ k ] =
mic_CSH [ i ] [ j ] [ k ];
12896 ArrPerc [ i ] [ j ] [ k ] = 0;
12902 for ( k = 0; k < 51; k++ ) {
12911 for ( k = 0; k <
SYSIZE; k++ ) {
12912 for ( j = 0; j <
SYSIZE; j++ ) {
12915 for ( kh = 0; kh < 51; kh++ ) {
12916 phase_temp [ kh ] = 0;
12921 while (
last != NULL ) {
12929 px =
cx(i, j, k, d1, d2, d3);
12930 py =
cy(i, j, k, d1, d2, d3);
12931 pz =
cz(i, j, k, d1, d2, d3);
12933 if (
IsSolidPhase(newmat [ px ] [ py ] [ pz ]) == 1 ) {
12935 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
12939 newmat [ px ] [ py ] [ pz ] =
BURNT;
12946 nmatx [ ncur ] = i;
12947 nmaty [ ncur ] = j;
12948 nmatz [ ncur ] = k;
12952 for ( inew = 1; inew <= ncur; inew++ ) {
12953 xcn = nmatx [ inew ];
12954 ycn = nmaty [ inew ];
12955 zcn = nmatz [ inew ];
12957 qx =
cx(xcn, ycn, zcn, d1, d2, d3);
12958 qy =
cy(xcn, ycn, zcn, d1, d2, d3);
12959 qz =
cz(xcn, ycn, zcn, d1, d2, d3);
12962 for ( jnew = 1; jnew <= 6; jnew++ ) {
12991 if ( y1 >= SYSIZE ) {
12993 }
else if ( y1 < 0 ) {
12998 if ( z1 >= SYSIZE ) {
13000 }
else if ( z1 < 0 ) {
13005 if ( ( x1 >= 0 ) && ( x1 < SYSIZE ) ) {
13006 px =
cx(x1, y1, z1, d1, d2, d3);
13007 py =
cy(x1, y1, z1, d1, d2, d3);
13008 pz =
cz(x1, y1, z1, d1, d2, d3);
13012 if ( (
IsSolidPhase(newmat [ px ] [ py ] [ pz ]) == 1 ) &&
13013 ( newmat [ px ] [ py ] [ pz ] !=
C3S ) &&
13014 ( newmat [ px ] [ py ] [ pz ] !=
C2S ) &&
13015 ( newmat [ px ] [ py ] [ pz ] !=
C3A ) &&
13016 ( newmat [ px ] [ py ] [ pz ] !=
C4AF ) ) {
13020 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
13022 newmat [ px ] [ py ] [ pz ] =
BURNT;
13025 printf(
"error in size of nnew %ld\n", nnew);
13028 nnewx [ nnew ] = x1;
13029 nnewy [ nnew ] = y1;
13030 nnewz [ nnew ] = z1;
13038 ( ( newmat [ px ] [ py ] [ pz ] ==
C3S ) ||
13039 ( newmat [ px ] [ py ] [ pz ] ==
C2S ) ||
13040 ( newmat [ px ] [ py ] [ pz ] ==
C3A ) ||
13041 ( newmat [ px ] [ py ] [ pz ] ==
C4AF ) ) ) {
13042 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
13044 newmat [ px ] [ py ] [ pz ] =
BURNT;
13047 printf(
"error in size of nnew %ld\n", nnew);
13050 nnewx [ nnew ] = x1;
13051 nnewy [ nnew ] = y1;
13052 nnewz [ nnew ] = z1;
13056 else if ( (
micpart [ qx ] [ qy ] [ qz ] ==
micpart [ px ] [ py ] [ pz ] ) &&
13057 ( ( newmat [ px ] [ py ] [ pz ] ==
C3S ) ||
13058 ( newmat [ px ] [ py ] [ pz ] ==
C2S ) ||
13059 ( newmat [ px ] [ py ] [ pz ] ==
C3A ) ||
13060 ( newmat [ px ] [ py ] [ pz ] ==
C4AF ) ) &&
13061 ( (
mic_CSH [ qx ] [ qy ] [ qz ] ==
C3S ) ||
13066 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
13068 newmat [ px ] [ py ] [ pz ] =
BURNT;
13071 printf(
"error in size of nnew %ld\n", nnew);
13074 nnewx [ nnew ] = x1;
13075 nnewy [ nnew ] = y1;
13076 nnewz [ nnew ] = z1;
13090 for ( icur = 1; icur <= ncur; icur++ ) {
13091 nmatx [ icur ] = nnewx [ icur ];
13092 nmaty [ icur ] = nnewy [ icur ];
13093 nmatz [ icur ] = nnewz [ icur ];
13096 }
while ( nnew > 0 );
13102 for ( j1 = 0; j1 <
SYSIZE; j1++ ) {
13103 for ( k1 = 0; k1 <
SYSIZE; k1++ ) {
13104 px =
cx(xl, j1, k1, d1, d2, d3);
13105 py =
cy(xl, j1, k1, d1, d2, d3);
13106 pz =
cz(xl, j1, k1, d1, d2, d3);
13107 qx =
cx(xh, j1, k1, d1, d2, d3);
13108 qy =
cy(xh, j1, k1, d1, d2, d3);
13109 qz =
cz(xh, j1, k1, d1, d2, d3);
13110 if ( ( newmat [ px ] [ py ] [ pz ] ==
BURNT ) && ( newmat [ qx ] [ qy ] [ qz ] ==
BURNT ) ) {
13114 if ( newmat [ px ] [ py ] [ pz ] ==
BURNT ) {
13115 newmat [ px ] [ py ] [ pz ] =
BURNT + 1;
13118 if ( newmat [ qx ] [ qy ] [ qz ] ==
BURNT ) {
13119 newmat [ qx ] [ qy ] [ qz ] =
BURNT + 1;
13124 if ( igood == 2 ) {
13127 while (
last != NULL ) {
13137 for ( kh = 0; kh < 51; kh++ ) {
13138 phase [ kh ] += phase_temp [ kh ];
13159 for ( kh = 0; kh < 51; kh++ ) {
13160 cnt_perc +=
phase [ kh ];
13170 for ( kh = 1; kh <=
HDCSH; kh++ ) {
13175 cnt_tot +=
count [ kh ];
13179 fprintf(
perc_phases,
"%d %f | ", cnt_tot, (
double ) cnt_perc / cnt_tot);
13180 for ( kh = 0; kh <=
HDCSH; kh++ ) {
13215 int CentPhase, NeighPhase;
13225 if ( ( NeighPhase !=
C3S ) &&
13226 NeighPhase !=
C2S &&
13227 NeighPhase !=
C3A &&
13228 NeighPhase !=
C4AF ) {
13233 else if ( ( CentPhase !=
C3S &&
13234 CentPhase !=
C2S &&
13235 CentPhase !=
C3A &&
13236 CentPhase !=
C4AF ) &&
13237 ( NeighPhase ==
C3S ||
13238 NeighPhase ==
C2S ||
13239 NeighPhase ==
C3A ||
13240 NeighPhase ==
C4AF ) ) {
13248 ( CentPhase ==
C3S ||
13249 CentPhase ==
C2S ||
13250 CentPhase ==
C3A ||
13251 CentPhase ==
C4AF ) &&
13252 ( NeighPhase ==
C3S ||
13253 NeighPhase ==
C2S ||
13254 NeighPhase ==
C3A ||
13255 NeighPhase ==
C4AF ) ) {
13311 for ( cz = 0; cz <
SYSIZE; cz++ ) {
13312 for ( cy = 0; cy <
SYSIZE; cy++ ) {
13313 for ( cx = 0; cx <
SYSIZE; cx++ ) {
13432 IsConnected(cx + 1, cy, cz + 1, 0, 1, 0) != 2 ) &&
13435 IsConnected(cx + 1, cy + 1, cz, 0, 0, 1) != 2 ) &&
13439 IsConnected(cx + 1, cy + 1, cz, 0, 0, 1) != 2 ) &&
13442 IsConnected(cx, cy + 1, cz + 1, 1, 0, 0) != 2 ) &&
13446 IsConnected(cx + 1, cy, cz + 1, 0, 1, 0) != 2 ) &&
13449 IsConnected(cx, cy + 1, cz + 1, 1, 0, 0) != 2 ) ) {
13463 IsConnected(cx + 1, cy, cz - 1, 0, 1, 0) != 2 ) &&
13466 IsConnected(cx + 1, cy + 1, cz, 0, 0, -1) != 2 ) &&
13470 IsConnected(cx + 1, cy + 1, cz, 0, 0, -1) != 2 ) &&
13473 IsConnected(cx, cy + 1, cz - 1, 1, 0, 0) != 2 ) &&
13477 IsConnected(cx + 1, cy, cz - 1, 0, 1, 0) != 2 ) &&
13480 IsConnected(cx, cy + 1, cz - 1, 1, 0, 0) != 2 ) ) {
13494 IsConnected(cx - 1, cy, cz - 1, 0, 1, 0) != 2 ) &&
13497 IsConnected(cx - 1, cy + 1, cz, 0, 0, -1) != 2 ) &&
13501 IsConnected(cx - 1, cy + 1, cz, 0, 0, -1) != 2 ) &&
13504 IsConnected(cx, cy + 1, cz - 1, -1, 0, 0) != 2 ) &&
13508 IsConnected(cx - 1, cy, cz - 1, 0, 1, 0) != 2 ) &&
13511 IsConnected(cx, cy + 1, cz - 1, -1, 0, 0) != 2 ) ) {
13524 IsConnected(cx - 1, cy, cz + 1, 0, 1, 0) != 2 ) &&
13527 IsConnected(cx - 1, cy + 1, cz, 0, 0, 1) != 2 ) &&
13531 IsConnected(cx - 1, cy + 1, cz, 0, 0, 1) != 2 ) &&
13534 IsConnected(cx, cy + 1, cz + 1, -1, 0, 0) != 2 ) &&
13538 IsConnected(cx - 1, cy, cz + 1, 0, 1, 0) != 2 ) &&
13541 IsConnected(cx, cy + 1, cz + 1, -1, 0, 0) != 2 ) ) {
13562 char extension [ 10 ];
13565 prefix = (
char * ) malloc(80);
13572 sprintf(extension,
"%04d",
icyc);
13573 strcpy(prefix,
"perc/out5.");
13574 strcat(prefix, extension);
13575 strcat(prefix,
".p.img");
13577 printf(
"Name of percolated output file is %s\n", prefix);
13581 if ( ( perc_img = fopen(prefix,
"w") ) == NULL ) {
13582 printf(
"\nFile %s can not be opened\n", prefix);
13587 for (
int dz = 0; dz <
SYSIZE; dz++ ) {
13588 for (
int dy = 0; dy <
SYSIZE; dy++ ) {
13589 for (
int dx = 0; dx <
SYSIZE; dx++ ) {
13591 fprintf(perc_img,
"%d\n",
ArrPerc [ dx ] [ dy ] [ dz ]);
13598 printf(
"Percolated file %s wrote\n", prefix);
13608 if (
last != NULL ) {
13629 if ( coord >=
SYSIZE ) {
13666 if ( * p_arr ==
CSH ) {
13690 CSH_vicinity [ i ] = 0;
13694 for ( cx = 0; cx <
SYSIZE; cx++ ) {
13695 for ( cy = 0; cy <
SYSIZE; cy++ ) {
13696 for ( cz = 0; cz <
SYSIZE; cz++ ) {
13704 for ( cx = 0; cx <
SYSIZE; cx++ ) {
13705 for ( cy = 0; cy <
SYSIZE; cy++ ) {
13706 for ( cz = 0; cz <
SYSIZE; cz++ ) {
13709 TotPhase +=
NumSol(cx, cy, cz);
13723 printf(
"\nNumerical Recipes run-time error...\n");
13724 printf(
"%s\n", error_text);
13725 printf(
"...now exiting to system...\n");
13734 v = (
float * ) malloc( (
unsigned ) ( nh - nl + 1 ) *
sizeof(
float ) );
13736 nrerror(
"allocation failure in vector()");
13746 v = (
int * ) malloc( (
unsigned ) ( nh - nl + 1 ) *
sizeof(
int ) );
13748 nrerror(
"allocation failure in ivector()");
13758 v = (
double * ) malloc( (
unsigned ) ( nh - nl + 1 ) *
sizeof(
double ) );
13760 nrerror(
"allocation failure in dvector()");
13774 m = (
float ** ) malloc( (
unsigned ) ( nrh - nrl + 1 ) *
sizeof(
float * ) );
13776 nrerror(
"allocation failure 1 in matrix_cem()");
13781 for ( i = nrl; i <= nrh; i++ ) {
13782 m [ i ] = (
float * ) malloc( (
unsigned ) ( nch - ncl + 1 ) *
sizeof(
float ) );
13784 nrerror(
"allocation failure 2 in matrix_cem()");
13798 m = (
double ** ) malloc( (
unsigned ) ( nrh - nrl + 1 ) *
sizeof(
double * ) );
13800 nrerror(
"allocation failure 1 in dmatrix()");
13805 for ( i = nrl; i <= nrh; i++ ) {
13806 m [ i ] = (
double * ) malloc( (
unsigned ) ( nch - ncl + 1 ) *
sizeof(
double ) );
13808 nrerror(
"allocation failure 2 in dmatrix()");
13821 m = (
int ** ) malloc( (
unsigned ) ( nrh - nrl + 1 ) *
sizeof(
int * ) );
13823 nrerror(
"allocation failure 1 in imatrix()");
13828 for ( i = nrl; i <= nrh; i++ ) {
13829 m [ i ] = (
int * ) malloc( (
unsigned ) ( nch - ncl + 1 ) *
sizeof(
int ) );
13831 nrerror(
"allocation failure 2 in imatrix()");
13847 m = (
float ** ) malloc( (
unsigned ) ( oldrh - oldrl + 1 ) *
sizeof(
float * ) );
13849 nrerror(
"allocation failure in submatrix()");
13854 for ( i = oldrl, j = newrl; i <= oldrh; i++, j++ ) {
13855 m [ j ] = a [ i ] + oldcl - newcl;
13873 for ( i = nrh; i >= nrl; i-- ) {
13874 free( (
char * ) ( m [ i ] + ncl ) );
13877 free( (
char * ) ( m + nrl ) );
13884 for ( i = nrh; i >= nrl; i-- ) {
13885 free( (
char * ) ( m [ i ] + ncl ) );
13888 free( (
char * ) ( m + nrl ) );
13895 for ( i = nrh; i >= nrl; i-- ) {
13896 free( (
char * ) ( m [ i ] + ncl ) );
13899 free( (
char * ) ( m + nrl ) );
13904 free( (
char * ) ( b + nrl ) );
13909 int i, j, nrow, ncol;
13912 nrow = nrh - nrl + 1;
13913 ncol = nch - ncl + 1;
13914 m = (
float ** ) malloc( (
unsigned ) ( nrow ) *
sizeof(
float * ) );
13916 nrerror(
"allocation failure in convert_matrix()");
13920 for ( i = 0, j = nrl; i <= nrow - 1; i++, j++ ) {
13921 m [ j ] = a + ncol * i - ncl;
13929 free( (
char * ) ( b + nrl ) );
13951 c.
r = a.
r * b.
r - a.
i * b.
i;
13952 c.
i = a.
i * b.
r + a.
r * b.
i;
13976 if ( fabs(b.
r) >= fabs(b.
i) ) {
13978 den = b.
r + r * b.
i;
13979 c.
r = ( a.
r + r * a.
i ) / den;
13980 c.
i = ( a.
i - r * a.
r ) / den;
13983 den = b.
i + r * b.
r;
13984 c.
r = ( a.
r * r + a.
i ) / den;
13985 c.
i = ( a.
i * r - a.
r ) / den;
13993 float x, y, ans, temp;
13998 }
else if ( y == 0.0 ) {
14000 }
else if ( x > y ) {
14002 ans = x * sqrt(1.0 + temp * temp);
14005 ans = y * sqrt(1.0 + temp * temp);
14015 if ( ( z.
r == 0.0 ) && ( z.
i == 0.0 ) ) {
14024 w = sqrt(x) * sqrt( 0.5 * ( 1.0 + sqrt(1.0 + r * r) ) );
14027 w = sqrt(y) * sqrt( 0.5 * ( r + sqrt(1.0 + r * r) ) );
14030 if ( z.
r >= 0.0 ) {
14032 c.
i = z.
i / ( 2.0 * w );
14034 c.
i = ( z.
i >= 0 ) ? w : -w;
14035 c.
r = z.
i / ( 2.0 * c.
i );
14075 sprintf(my_string,
"%d",
iseed);
14080 return ( 1 - b - c ) * x + ( 1 - a - c ) * y + ( 1 - a - b ) * z;
14086 return ( 1 - a - b ) * x + ( 1 - b - c ) * y + ( 1 - a - c ) * z;
14091 return ( 1 - a - c ) * x + ( 1 - a - b ) * y + ( 1 - b - c ) * z;
14102 double E_CSH_hmg, nu_CSH_hmg;
14108 for ( x = 0; x < 34; x++ ) {
14112 for ( x = 0; x < 51; x++ ) {
14114 if ( x ==
HDCSH ) {
14116 }
else if ( x ==
EMPTYP ) {
14123 if ( perc_unperc_flag == 0 ) {
14134 for ( x = 2; x < 34; x++ ) {
14143 for ( x = 0; x < 34; x++ ) {
14149 PhaseMatrix(0, 1) = 0.001;
14150 PhaseMatrix(0, 2) = 0.499924;
14152 PhaseMatrix(1, 1) = 135.0;
14153 PhaseMatrix(1, 2) = 0.300000;
14155 PhaseMatrix(2, 1) = 130.0;
14156 PhaseMatrix(2, 2) = 0.300000;
14158 PhaseMatrix(3, 1) = 145.0;
14159 PhaseMatrix(3, 2) = 0.300000;
14161 PhaseMatrix(4, 1) = 125.0;
14162 PhaseMatrix(4, 2) = 0.300000;
14164 PhaseMatrix(5, 1) = 30.00;
14165 PhaseMatrix(5, 2) = 0.300000;
14167 PhaseMatrix(6, 1) = 62.90;
14168 PhaseMatrix(6, 2) = 0.300000;
14170 PhaseMatrix(7, 1) = 73.60;
14171 PhaseMatrix(7, 2) = 0.275000;
14173 PhaseMatrix(8, 1) = 72.80;
14174 PhaseMatrix(8, 2) = 0.167000;
14176 PhaseMatrix(9, 1) = 79.60;
14177 PhaseMatrix(9, 2) = 0.309900;
14179 PhaseMatrix(10, 1) = 72.80;
14180 PhaseMatrix(10, 2) = 0.16700;
14182 PhaseMatrix(11, 1) = 72.80;
14183 PhaseMatrix(11, 2) = 0.16700;
14185 PhaseMatrix(12, 1) = 72.80;
14186 PhaseMatrix(12, 2) = 0.16700;
14188 PhaseMatrix(13, 1) = 38.00;
14189 PhaseMatrix(13, 2) = 0.30500;
14190 PhaseMatrix(14, 0) = -10.;
14191 PhaseMatrix(14, 1) = -10.0;
14192 PhaseMatrix(14, 2) = -10.000;
14194 PhaseMatrix(15, 1) = 22.40;
14195 PhaseMatrix(15, 2) = 0.25000;
14197 PhaseMatrix(16, 1) = 22.40;
14198 PhaseMatrix(16, 2) = 0.25000;
14200 PhaseMatrix(17, 1) = 22.40;
14201 PhaseMatrix(17, 2) = 0.25000;
14203 PhaseMatrix(18, 1) = 42.30;
14204 PhaseMatrix(18, 2) = 0.32380;
14206 PhaseMatrix(19, 1) = 22.40;
14207 PhaseMatrix(19, 2) = 0.24590;
14209 PhaseMatrix(20, 1) = 22.40;
14210 PhaseMatrix(20, 2) = 0.24590;
14212 PhaseMatrix(21, 1) = 22.40;
14213 PhaseMatrix(21, 2) = 0.24590;
14215 PhaseMatrix(22, 1) = 42.30;
14216 PhaseMatrix(22, 2) = 0.32380;
14218 PhaseMatrix(23, 1) = 22.40;
14219 PhaseMatrix(23, 2) = 0.24590;
14221 PhaseMatrix(24, 1) = 22.40;
14222 PhaseMatrix(24, 2) = 0.25000;
14224 PhaseMatrix(25, 1) = 30.00;
14225 PhaseMatrix(25, 2) = 0.30000;
14227 PhaseMatrix(26, 1) = 60.00;
14228 PhaseMatrix(26, 2) = 0.30000;
14230 PhaseMatrix(27, 1) = 42.30;
14231 PhaseMatrix(27, 2) = 0.32380;
14233 PhaseMatrix(28, 1) = 79.60;
14234 PhaseMatrix(28, 2) = 0.30990;
14236 PhaseMatrix(29, 1) = 30.00;
14237 PhaseMatrix(29, 2) = 0.30000;
14239 PhaseMatrix(30, 1) = 0.001;
14240 PhaseMatrix(30, 2) = 0.00100;
14251 LevelI(0, 1) = 21.7;
14252 LevelI(0, 2) = 0.24;
14254 LevelI(1, 1) = 29.4;
14255 LevelI(1, 2) = 0.24;
14261 for ( i = 0; i < 31; i++ ) {
14262 PhaseMatrix(i, 0) =
PhaseFrac [ i + 1 ];
14268 PhaseMatrix(14, 1) = CSH_level.
E_hmg;
14269 PhaseMatrix(14, 2) = CSH_level.
nu_hmg;
14271 E_CSH_hmg = CSH_level.
E_hmg;
14272 nu_CSH_hmg = CSH_level.
nu_hmg;
14276 E = Paste_level.
E_hmg;
14277 nu = Paste_level.
nu_hmg;
14279 ( void ) E_CSH_hmg;
14280 ( void ) nu_CSH_hmg;
14286 Homogenize Paste_level, Mortar_level, Concrete_level;
14290 FloatMatrix PhaseMatrix(4, 3), Mortar(4, 3), Concrete(4, 3);
14293 PhaseMatrix(0, 1) = E_paste_inp;
14294 PhaseMatrix(0, 2) = nu_paste_inp;
14301 PhaseMatrix(3, 0) = Vol_entrained_entrapped_air / vol_tot_paste;
14302 PhaseMatrix(3, 1) = 0.001;
14303 PhaseMatrix(3, 2) = 0.001;
14306 * E_paste = Paste_level.
E_hmg;
14307 * nu_paste = Paste_level.
nu_hmg;
14312 double n_FA, n_CA, vol_ITZ_FA, vol_ITZ_CA;
14317 double vol_tot_mortar = vol_tot_paste +
Vol_FA;
14319 Mortar(0, 0) = Vol_FA / vol_tot_mortar;
14322 Mortar(1, 0) = vol_ITZ_FA / vol_tot_mortar;
14324 Mortar(1, 2) = * nu_paste;
14325 Mortar(2, 0) = ( vol_tot_paste - vol_ITZ_FA ) / vol_tot_mortar;
14326 Mortar(2, 1) = * E_paste;
14327 Mortar(2, 2) = * nu_paste;
14329 Mortar(3, 1) = 10.;
14330 Mortar(3, 2) = 0.3;
14333 * E_mortar = Mortar_level.
E_hmg;
14334 * nu_mortar = Mortar_level.
nu_hmg;
14337 double vol_tot_concrete = vol_tot_mortar +
Vol_CA;
14339 Concrete(0, 0) = Vol_CA / vol_tot_concrete;
14342 Concrete(1, 0) = vol_ITZ_CA / vol_tot_concrete;
14344 Concrete(1, 2) = * nu_paste;
14345 Concrete(2, 0) = ( vol_tot_mortar - vol_ITZ_CA ) / vol_tot_concrete;
14346 Concrete(2, 1) = * E_mortar;
14347 Concrete(2, 2) = * nu_mortar;
14348 Concrete(3, 0) = 0.;
14349 Concrete(3, 1) = 10.;
14350 Concrete(3, 2) = 0.3;
14353 E_concrete = Concrete_level.
E_hmg;
14354 nu_concrete = Concrete_level.
nu_hmg;
14363 c3s = ( double )
c3sinit / sum;
14364 c2s = ( double )
c2sinit / sum;
14365 c3a = ( double )
c3ainit / sum;
14367 gypsum = ( double )
ncsbar / sum;
14368 hemi = ( double )
heminit / sum;
14369 anh = ( double ) anhinit / sum;
14374 #ifdef __TM_MODULE //transport module for OOFEM 14393 if ( !cemhydmat->
eachGP ) {
14413 fprintf(file,
" status {");
14414 fprintf( file,
"cyc %d time %e DoH %f conductivity %f capacity %f density %f", cemhydmat->
giveCycleNumber(this->gp), 3600 * cemhydmat->
giveTimeOfCycle(this->
gp), cemhydmat->
giveDoHActual(this->
gp), cemhydmat->
giveIsotropicConductivity(this->
gp,tStep), cemhydmat->
giveConcreteCapacity(this->
gp, tStep), cemhydmat->
giveConcreteDensity(this->
gp, tStep) );
14420 if ( !cemhydmat->
eachGP ) {
14422 fprintf( file,
" master material %d", cemhydmat->
giveNumber() );
14424 fprintf( file,
" slave of material %d", cemhydmat->
giveNumber() );
14427 fprintf( file,
" independent microstructure %p from material %d",
this, cemhydmat->
giveNumber() );
14430 fprintf(file,
"}\n");
14442 printf(
"Constructor of CemhydMatStatus called\n");
14449 CemhydMatStatus :: InitializePy(
const char *inp)
14461 #include <boost/python.hpp> 14462 BOOST_PYTHON_MODULE(cemhydmodule)
14465 boost :: python :: class_< oofem :: CemhydMatStatus >(
"cemhydmatstatus")
14466 .def(
"InitializePy", & oofem :: CemhydMatStatus :: InitializePy)
int movehem(int xcur, int ycur, int zcur, int finalstep, float nucprgyp)
void free_vector(float *v, int nl)
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
int moveas(int xcur, int ycur, int zcur, int finalstep)
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
int burnset(int d1, int d2, int d3)
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in full form.
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
void outputImageFileUnperc(char ***m)
void GetInputParams(char *my_string)
int movecaco3(int xcur, int ycur, int zcur, int finalstep)
int burn3d(int npix, int d1, int d2, int d3)
double giveTotalVolume(void)
double GiveDoHLastCyc(void)
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
Access method for default integration rule.
int gsphere(int numgen, long int *numeach, int *sizeeach, int *pheach)
void dealloc_int_3D(int ***(&mask), long SYSIZE)
int ** imatrix(int nrl, int nrh, int ncl, int nch)
double Vol_cement_clinker_gypsum
fcomplex_cem Conjg(fcomplex_cem z)
int IsSolidPhase(int phase)
#define _IFT_CemhydMat_densitytype
int readInputFileAndInitialize(const char *inp, bool generateMicrostructure)
void free_dmatrix(double **m, int nrl, int nrh, int ncl)
Domain * domain
Link to domain object, useful for communicating with other FEM components.
int moveone(int *xloc, int *yloc, int *zloc, int *act, int sumold)
std::string XMLfileName
XML input file name for CEMHYD3D.
void addrand(int randid, long int nneed)
unsigned int * CSH_vicinity
void free_imatrix(int **m, int nrl, int nrh, int ncl)
Class for elastic homogenization.
int extettr(int xpres, int ypres, int zpres, int etype)
void nrerror(const char *error_text)
void zero()
Sets all component to zero.
double & at(int i)
Coefficient access function.
struct percolatedpath * last
double PartHeat
The last incremental heat returned from a GP.
void extpozz(int xpres, int ypres, int zpres)
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
CemhydMatStatus(int n, Domain *d, GaussPoint *gp, CemhydMatStatus *CemStat, CemhydMat *cemhydmat, bool withMicrostructure)
Create status in an integration point.
int eachGP
Assign a separate microstructure in each integration point.
int chksph(int xin, int yin, int zin, int radd, int wflg, int phasein, int phase2)
fcomplex_cem Csqrt(fcomplex_cem z)
void hydrate(int fincyc, int stepmax, float chpar1, float chpar2, float hgpar1, float hgpar2, float fhpar1, float fhpar2, float gypar1, float gypar2)
void free_ivector(int *v, int nl)
int MoveToTime(double GiveTemp, double TargTime)
#define _IFT_CemhydMat_conductivitytype
void GetInitClinkerPhases(double &c3s, double &c2s, double &c3a, double &c4af, double &gypsum, double &hemi, double &anh)
virtual void updateYourself(TimeStep *tStep)
Update equilibrium history variables according to temp-variables.
int movepix(int ntomove, int ph1, int ph2)
void AnalyticHomogenizationPaste(double &E, double &nu, int perc_unperc_flag)
const char * __MatResponseModeToString(MatResponseMode _value)
int movecacl2(int xcur, int ycur, int zcur, int finalstep)
int * ivector(int nl, int nh)
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver's output to given stream.
int movecsh(int xcur, int ycur, int zcur, int finalstep, int cycorig)
double giveTargetTime()
Returns target time.
void free_convert_matrix(float **b, int nrl)
void WriteUnsortedList(int px, int py, int pz)
Abstract base class for all finite elements.
double giveAverageTemperature(void)
double Mass_cement_concrete
Element * giveElement()
Returns corresponding element to receiver.
double averageTemperature
Average temperature through integration points.
This class implements a transport material status information.
int countbox(int boxsize, int qx, int qy, int qz)
long cy(int x, int y, int z, int a, int b, int c)
void alloc_char_3D(char ***(&mic), long SYSIZE)
int IsConnected(int cx, int cy, int cz, int dx, int dy, int dz)
fcomplex_cem Cadd(fcomplex_cem a, fcomplex_cem b)
double * dvector(int nl, int nh)
int & at(int i)
Coefficient access function.
MatResponseMode
Describes the character of characteristic material matrix.
struct percolatedpath * prev
int reinforcementDegree
Degree of reinforcement, if defined, reinforcement effect for conductivity and capacity is accounted ...
void initializeMicrostructure(void)
int movegyp(int xcur, int ycur, int zcur, int finalstep)
void alloc_shortint_3D(short int ***(&mic), long SYSIZE)
void free_submatrix(float *b, int nrl)
int countboxc(int boxsize, int qx, int qy, int qz)
double MoveCycles(double GiveTemp, int cycles)
virtual double giveDoHActual(GaussPoint *gp)
Returns DoH of the closest CEMHYD3D cycle after the target time.
FloatArray scaling
Array containing scaling factors for density, conductivity and capacity.
double Concrete_bulk_density
double init_material_time
Inital material time for growing problems.
float ** convert_matrix(float *a, int nrl, int nrh, int ncl, int nch)
virtual double giveConcreteCapacity(GaussPoint *gp, TimeStep *tStep)
Returns concrete thermal capacity depending on chosen type.
void extafm(int xpres, int ypres, int zpres)
void zroots(fcomplex_cem a[], int m, fcomplex_cem roots[], int polish)
void QueryNumAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, int &val)
virtual double computeVolumeAround(GaussPoint *gp)
Returns volume related to given integration point.
IntArray nowarnings
Array containing warnings supression for density, conductivity, capacity, high temperature.
void extfh3(int xpres, int ypres, int zpres)
int chkfloc(int xin, int yin, int zin, int radd)
int surfpix(int xin, int yin, int zin)
GaussPoint * gp
Stores GP of the CemhydMatStatus.
void disrealnew_init(void)
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
CemhydMatStatus * MasterCemhydMatStatus
Pointer to master CemhydMatStatus, which is shared among related integration points (on element...
const FloatArray & giveField()
Return last field.
void sysscan(int ph1, int ph2)
int conductivityType
Use different methods to evaluate material parameters.
long cz(int x, int y, int z, int a, int b, int c)
struct percolatedpath * current
virtual void averageTemperature()
Perform averaging on a master CemhydMatStatus.
int movec4a(int xcur, int ycur, int zcur, int finalstep, float nucprob)
void dealloc_long_3D(long ***(&mic), long SYSIZE)
void GenerateConnNumbers(void)
#define _IFT_CemhydMat_nowarnings
void alloc_long_3D(long ***(&mic), long SYSIZE)
int movec3a(int xcur, int ycur, int zcur, int finalstep, float nucprob)
int edgecnt(int xck, int yck, int zck, int ph1, int ph2, int ph3)
void outputImageFilePerc(void)
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
This class implements an isotropic linear heat material.
virtual double giveIsotropicConductivity(GaussPoint *gp, TimeStep *tStep)
Returns concrete heat conductivity depending on chosen type.
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Print receiver's output to given stream.
Material * giveMaterial()
Returns reference to material associated to related element of receiver.
void dealloc_double_3D(double ***(&mic), long SYSIZE)
void extc3ah6(int xpres, int ypres, int zpres)
int extstrat(int xpres, int ypres, int zpres)
void resize(int n)
Checks size of receiver towards requested bounds.
int extfreidel(int xpres, int ypres, int zpres)
fcomplex_cem RCmul(float x, fcomplex_cem a)
virtual ~CemhydMat()
Destructor.
virtual double giveConcreteDensity(GaussPoint *gp, TimeStep *tStep)
Returns concrete density depending on chosen type.
Abstract base class representing a material status information.
void dealloc_char_3D(char ***(&mic), long SYSIZE)
int movecas2(int xcur, int ycur, int zcur, int finalstep)
int moveettr(int xcur, int ycur, int zcur, int finalstep)
Class representing vector of real numbers.
void selfConsistent(FloatMatrix &PhaseMatrix)
Self-consistent homogenization method of Hill and Budiansky for spherical isotropic inclusions...
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
int Calculate_elastic_homogenization
Flag to proceed percolation filtering and elastic homogenization.
Implementation of matrix containing floating point numbers.
double nu_hmg
Effective Poisson's ratio.
IRResultType
Type defining the return values of InputRecord reading operations.
void burn_phases(int d1, int d2, int d3)
double E_hmg
Effective Young's modulus or the lower bound.
void makeinert(long int ndesire)
int loccsh(int xcur, int ycur, int zcur, int extent)
struct percolatedpath * next
void free_dvector(double *v, int nl)
#define _IFT_CemhydMat_capacitytype
void laguer(fcomplex_cem a[], int m, fcomplex_cem *x, float eps, int polish)
void moriTanaka(FloatMatrix &PhaseMatrix, int refRow)
Mori-Tanaka homogenization method for spherical isotropic inclusions.
virtual int initMaterial(Element *element)
Optional function to call specific procedures when initializing a material.
void disrealnew(double GiveTemp, double TargTime, int flag)
void AnalyticHomogenizationConcrete(double E_paste_inp, double nu_paste_inp, double *E_paste, double *nu_paste, double *E_mortar, double *nu_mortar, double &E_concrete, double &nu_concrete)
void sysinit(int ph1, int ph2)
float Cabs(fcomplex_cem z)
#define _IFT_CemhydMat_inputFileName
void PercolateForOutput(void)
double ** dmatrix(int nrl, int nrh, int ncl, int nch)
void QueryStringAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, char *chars)
void alloc_double_3D(double ***(&mic), long SYSIZE)
long int target_anhydrite
double Concrete_thermal_conductivity
float ** submatrix(float **a, int oldrl, int oldrh, int oldcl, int oldch, int newrl, int newcl)
void rand3d(int phasein, int phaseout, float xpt)
int MoveToDoH(double GiveTemp, double DesiredDoH, int maxcyc)
double GiveDoHActual(void)
Return degree of hydration of the receiver.
virtual double give(int aProperty, GaussPoint *gp, TimeStep *tStep)
int countem(int xp, int yp, int zp, int phin)
void dealloc_shortint_3D(short int ***(&mic), long SYSIZE)
int moveanh(int xcur, int ycur, int zcur, int finalstep, float nucprgyp)
void passone(int low, int high, int cycid, int cshexflag)
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual double giveTimeOfCycle(GaussPoint *gp)
Returns time of the CEMHYD3D at the first cycle after the target time.
virtual void storeWeightTemperatureProductVolume(Element *element, TimeStep *tStep)
Store temperatures multiplied with volume around GPs - need before temperature averaging.
virtual ~CemhydMatStatus()
void sinter3d(int ph1id, int ph2id, float rhtarget)
int chckedge(int xck, int yck, int zck)
void extslagcsh(int xpres, int ypres, int zpres)
fcomplex_cem Csub(fcomplex_cem a, fcomplex_cem b)
fcomplex_cem Cmul(fcomplex_cem a, fcomplex_cem b)
CemhydMat(int n, Domain *d)
Constructor.
int movefh3(int xcur, int ycur, int zcur, int finalstep, float nucprob)
#define _IFT_CemhydMat_scaling
void drawfloc(int xin, int yin, int zin, int radd, int phasein, int phase2)
virtual void computeInternalSourceVector(FloatArray &val, GaussPoint *gp, TimeStep *tStep, ValueModeType mode)
Computes the internal source vector of receiver.
REGISTER_Material(DummyMaterial)
double castingTime
Casting time.
struct cluster * nextpart
float * vector(int nl, int nh)
bool isIcApply()
Check if receiver is solution step when initial conditions should apply.
virtual void clearWeightTemperatureProductVolume(Element *element)
Clear temperatures multiplied with volume around GPs - need before temperature averaging.
int giveSize() const
Returns the size of receiver.
long cx(int x, int y, int z, int a, int b, int c)
double GivePower(double GiveTemp, double TargTime)
the oofem namespace is to define a context or scope in which all oofem names are defined.
void setAverageTemperatureVolume(double temperature, double volume)
Auxiliary function for temperature averaging over GPs.
fcomplex_cem Cdiv(fcomplex_cem a, fcomplex_cem b)
double * last_values
Array for storing temporary values (elastic properties etc.)
float ** matrix_cem(int nrl, int nrh, int ncl, int nch)
int movech(int xcur, int ycur, int zcur, int finalstep, float nucprob)
void free_matrix(float **m, int nrl, int nrh, int ncl)
double computeConcreteCapacityBentz(void)
#define _IFT_CemhydMat_eachgp
void CSHbox(unsigned int *CSH_vicinity)
void alloc_int_3D(int ***(&mask), long SYSIZE)
Class representing integration point in finite element program.
#define OOFEM_WARNING(...)
double GiveTotCemHeat(void)
fcomplex_cem ComplexCemhyd(float re, float im)
int NumSol(int cx, int cy, int cz)
Class representing solution step.
double IPVolume
Volume associated to master IP of one CemhydMat.
void herveZaoui(FloatMatrix &PhaseMatrix)
Herve and Zaoui's homogenization scheme for n-spherical isotropic domains and arbitrary number of iso...
void extgyps(int xpres, int ypres, int zpres)
virtual double giveCharacteristicValue(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Compute heat thermal capacity per volume.
double Vol_entrained_entrapped_air
#define _IFT_CemhydMat_reinforcementDegree
const char * __ValueModeTypeToString(ValueModeType _value)
virtual int giveCycleNumber(GaussPoint *gp)
Returns cycle number at the closest cycle after the target time.
void resize(int s)
Resizes receiver towards requested size.
int icyc
Cycle of celular automata.