111 if ( mode == VM_Total || mode == VM_TotalIntrinsic ) {
119 OOFEM_ERROR(
"Undefined mode %s\n", __ValueModeTypeToString(mode) );
122 if ( mode == VM_Total || mode == VM_TotalIntrinsic ) {
125 OOFEM_ERROR(
"Undefined mode %s\n", __ValueModeTypeToString(mode) );
169 double conduct = 0.0;
176 conduct = IsotropicHeatTransferMaterial :: give(
'k', gp, tStep);
178 conduct = IsotropicHeatTransferMaterial :: give(
'k', gp, tStep) * ( 1.33 - 0.33 * ms->GiveDoHActual() );
187 if ( !this->
nowarnings.at(2) && ( conduct < 0.3 || conduct > 5 ) ) {
188 OOFEM_WARNING(
"Weird concrete thermal conductivity %f W/m/K\n", conduct);
191 conduct *= this->
scaling.at(2);
200 double capacityConcrete = 0.0;
207 capacityConcrete = IsotropicHeatTransferMaterial :: give(
'c', gp, tStep);
209 capacityConcrete = ms->computeConcreteCapacityBentz();
211 capacityConcrete = 1000 * ms->GiveCp();
219 if ( !this->
nowarnings.at(3) && ( capacityConcrete < 500 || capacityConcrete > 2000 ) ) {
220 OOFEM_WARNING(
"Weird concrete heat capacity %f J/kg/K\n", capacityConcrete);
223 capacityConcrete *= this->
scaling.at(3);
225 return capacityConcrete;
231 double concreteBulkDensity = 0.0;
238 concreteBulkDensity = IsotropicHeatTransferMaterial :: give(
'd', gp, tStep);
240 concreteBulkDensity = ms->GiveDensity();
248 if ( !this->
nowarnings.at(1) && ( concreteBulkDensity < 1000 || concreteBulkDensity > 4000 ) ) {
249 OOFEM_WARNING(
"Weird concrete density %f kg/m3\n", concreteBulkDensity);
252 concreteBulkDensity *= this->
scaling.at(1);
254 return concreteBulkDensity;
260 if ( mode == Capacity ) {
262 }
else if ( mode == IntSource ) {
266 double lastEquilibratedTemperature = status->
giveField();
268 double krate, EaOverR, val;
274 EaOverR = 1000. * ms->
E_act / 8.314;
276 if ( ms->
icyc > 1 ) {
277 krate = exp( -EaOverR * ( 1. / ( ms->
temp_cur + 273.15 ) - 1. / ( lastEquilibratedTemperature + 273.15 ) ) );
289 val = EaOverR * krate / ( ms->
temp_cur + 273.15 ) / ( ms->
temp_cur + 273.15 );
293 OOFEM_ERROR(
"unknown mode (%s)\n", __MatResponseModeToString(mode) );
310 if ( type == IST_HydrationDegree ) {
314 }
else if ( type == IST_Density ) {
318 }
else if ( type == IST_ThermalConductivityIsotropic ) {
322 }
else if ( type == IST_HeatCapacity ) {
326 }
else if ( type == IST_AverageTemperature ) {
330 }
else if ( type == IST_YoungModulusVirginPaste ) {
334 }
else if ( type == IST_PoissonRatioVirginPaste ) {
338 }
else if ( type == IST_YoungModulusConcrete ) {
342 }
else if ( type == IST_PoissonRatioConcrete ) {
347 return TransportMaterial :: giveIPValue(answer, gp, type, tStep);
366 gp->setMaterialStatus( ms );
373void CemhydMat :: clearWeightTemperatureProductVolume(
Element *element)
383void CemhydMat :: storeWeightTemperatureProductVolume(
Element *element,
TimeStep *tStep)
391 element->
giveIPValue(vecTemperature, gp, IST_Temperature, tStep);
397void CemhydMat :: averageTemperature()
409 IsotropicHeatTransferMaterial :: initializeFrom(ir);
418 for (
int i = 1; i <=
scaling.giveSize(); i++ ) {
434 if (
scaling.giveSize() != 3 ) {
446 OOFEM_ERROR(
"Use function CemhydMat :: initMaterial instead");
533 printf(
"Constructor of CemhydMatStatus on GP %p, withMicrostructure %d, copy from CemhydMatStatus %p\n",
gp, withMicrostructure, CemStat);
537 if ( withMicrostructure ) {
543 for ( k = 0; k <
SYSIZE; k++ ) {
544 for ( j = 0; j <
SYSIZE; j++ ) {
545 for ( i = 0; i <
SYSIZE; i++ ) {
548 mic [ i ] [ j ] [ k ] =
micorig [ i ] [ j ] [ k ];
558void CemhydMatStatus :: initializeMicrostructure()
822 phase =
new long int [ 51 ];
844 for (
int i = 0; i <=
EMPTYP; i++ ) {
849 for (
int i = 0; i <=
HDCSH; i++ ) {
860CemhydMatStatus :: ~CemhydMatStatus()
946 while (
headant->nextant != NULL ) {
951 printf(
"Deallocating headant\n");
961 printf(
"Deallocating arrays\n");
976void CemhydMatStatus :: alloc_char_3D(
char ***( &
mic ),
long SYSIZE)
979 for (
int x = 0; x <
SYSIZE; x++ ) {
981 for (
int y = 0; y <
SYSIZE; y++ ) {
983 if (
mic [ x ] [ y ] == NULL ) {
984 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
990void CemhydMatStatus :: dealloc_char_3D(
char ***( &
mic ),
long SYSIZE)
993 for (
int x = 0; x <
SYSIZE; x++ ) {
994 for (
int y = 0; y <
SYSIZE; y++ ) {
995 delete []
mic [ x ] [ y ];
1005void CemhydMatStatus :: alloc_long_3D(
long ***( &
mic ),
long SYSIZE)
1008 for (
int x = 0; x <
SYSIZE; x++ ) {
1010 for (
int y = 0; y <
SYSIZE; y++ ) {
1012 if (
mic [ x ] [ y ] == NULL ) {
1013 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1020void CemhydMatStatus :: dealloc_long_3D(
long ***( &
mic ),
long SYSIZE)
1022 if (
mic != NULL ) {
1023 for (
int x = 0; x <
SYSIZE; x++ ) {
1024 for (
int y = 0; y <
SYSIZE; y++ ) {
1025 delete []
mic [ x ] [ y ];
1028 delete []
mic [ x ];
1035void CemhydMatStatus :: alloc_int_3D(
int ***( &
mic ),
long SYSIZE)
1038 for (
int x = 0; x <
SYSIZE; x++ ) {
1040 for (
int y = 0; y <
SYSIZE; y++ ) {
1042 if (
mic [ x ] [ y ] == NULL ) {
1043 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1050void CemhydMatStatus :: dealloc_int_3D(
int ***( &
mic ),
long SYSIZE)
1052 if (
mic != NULL ) {
1053 for (
int x = 0; x <
SYSIZE; x++ ) {
1054 for (
int y = 0; y <
SYSIZE; y++ ) {
1055 delete []
mic [ x ] [ y ];
1058 delete []
mic [ x ];
1065void CemhydMatStatus :: alloc_shortint_3D(
short int ***( &
mic ),
long SYSIZE)
1068 for (
int x = 0; x <
SYSIZE; x++ ) {
1070 for (
int y = 0; y <
SYSIZE; y++ ) {
1071 mic [ x ] [ y ] =
new short int [
SYSIZE ];
1072 if (
mic [ x ] [ y ] == NULL ) {
1073 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1080void CemhydMatStatus :: dealloc_shortint_3D(
short int ***( &
mic ),
long SYSIZE)
1082 if (
mic != NULL ) {
1083 for (
int x = 0; x <
SYSIZE; x++ ) {
1084 for (
int y = 0; y <
SYSIZE; y++ ) {
1085 delete []
mic [ x ] [ y ];
1088 delete []
mic [ x ];
1095void CemhydMatStatus :: alloc_double_3D(
double ***( &
mic ),
long SYSIZE)
1098 for (
int x = 0; x <
SYSIZE; x++ ) {
1100 for (
int y = 0; y <
SYSIZE; y++ ) {
1101 mic [ x ] [ y ] =
new double [
SYSIZE ];
1102 if (
mic [ x ] [ y ] == NULL ) {
1103 printf(
"Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1110void CemhydMatStatus :: dealloc_double_3D(
double ***( &
mic ),
long SYSIZE)
1112 if (
mic != NULL ) {
1113 for (
int x = 0; x <
SYSIZE; x++ ) {
1114 for (
int y = 0; y <
SYSIZE; y++ ) {
1115 delete []
mic [ x ] [ y ];
1118 delete []
mic [ x ];
1127void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *
xmlFile,
const char *elementName,
int position,
int &val)
1131 XMLHandle docHandle = XMLHandle(
xmlFile);
1132 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1133 if ( elemSelected == NULL ) {
1134 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1138 sprintf(key,
"key%d", position);
1139 success = elemSelected->QueryIntAttribute(key, & val);
1140 if ( success != XML_SUCCESS ) {
1141 printf(
"Cannot read int value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1146void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *
xmlFile,
const char *elementName,
int position,
long int &val)
1150 val =
static_cast< long int >(temp);
1153void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *
xmlFile,
const char *elementName,
const char *key,
int &val)
1156 XMLHandle docHandle = XMLHandle(
xmlFile);
1157 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1158 if ( elemSelected == NULL ) {
1159 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1163 success = elemSelected->QueryIntAttribute(key, & val);
1164 if ( success != XML_SUCCESS ) {
1165 printf(
"Cannot read int value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1171void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *
xmlFile,
const char *elementName,
int position,
double &val)
1175 XMLHandle docHandle = XMLHandle(
xmlFile);
1176 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1177 if ( elemSelected == NULL ) {
1178 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1182 sprintf(key,
"key%d", position);
1183 success = elemSelected->QueryDoubleAttribute(key, & val);
1184 if ( success != XML_SUCCESS ) {
1185 printf(
"Cannot read double value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1190void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *
xmlFile,
const char *elementName,
const char *key,
double &val)
1193 XMLHandle docHandle = XMLHandle(
xmlFile);
1194 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1195 if ( elemSelected == NULL ) {
1196 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1200 success = elemSelected->QueryDoubleAttribute(key, & val);
1201 if ( success != XML_SUCCESS ) {
1202 printf(
"Cannot read double value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1207void CemhydMatStatus :: QueryStringAttributeExt(XMLDocument *
xmlFile,
const char *elementName,
int position,
char *chars)
1211 XMLHandle docHandle = XMLHandle(
xmlFile);
1213 XMLElement *elemSelected = docHandle.FirstChildElement(
"cemhyd").FirstChildElement(elementName).ToElement();
1214 if ( elemSelected == NULL ) {
1215 printf(
"Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1219 sprintf(key,
"key%d", position);
1223 const char *cstr = elemSelected->Attribute(key);
1225 success = XML_SUCCESS;
1227 success = XML_NO_ATTRIBUTE;
1229 if ( success != XML_SUCCESS ) {
1230 printf(
"Cannot read string value or key %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1233 str1 = std :: string(cstr);
1235 strcpy( chars, str1.c_str() );
1243int CemhydMatStatus :: readInputFileAndInitialize(
const char *inp,
bool generateMicrostructure)
1247 F =
new cmlfile(inp);
1250 F->set_label_string(0,
"Rand_seed_num");
1251 F->set_label_string(1,
"Input_img_file");
1252 F->set_label_string(2,
"Input_id_file");
1253 F->set_label_string(3,
"Saturated_sealed");
1254 F->set_label_string(4,
"Induction_time");
1255 F->set_label_string(5,
"Ea_cement");
1256 F->set_label_string(6,
"Ea_pozz");
1257 F->set_label_string(7,
"Ea_slag");
1258 F->set_label_string(8,
"Beta");
1259 F->set_label_string(9,
"Mass_SCM_FA_CA_inert_frac");
1260 F->set_label_string(10,
"Mass_cem");
1261 F->set_label_string(11,
"Cp_SCM_FA_CA_inert");
1262 F->set_label_string(12,
"Cp_cem");
1263 F->set_label_string(13,
"Given_microstructure");
1264 F->set_label_string(14,
"Output_initial_microstructure");
1265 F->set_label_string(15,
"Output_initial_microstructure_img_file");
1266 F->set_label_string(16,
"Output_initial_microstructure_id_file");
1267 F->set_label_string(17,
"Cycle_freq_perc_pore");
1268 F->set_label_string(18,
"Cycle_freq_perc_sol");
1269 F->set_label_string(19,
"Total_sodium");
1270 F->set_label_string(20,
"Total_potassium");
1271 F->set_label_string(21,
"Readily_soluble_sodium");
1272 F->set_label_string(22,
"Readily_soluble_potassium");
1273 F->set_label_string(23,
"Diffusion_steps_per_cycle");
1274 F->set_label_string(24,
"CH_nucleation_probability");
1275 F->set_label_string(25,
"CH_scale_factor");
1276 F->set_label_string(26,
"Gypsum_nucleation_probability");
1277 F->set_label_string(27,
"Gypsum_scale_factor");
1278 F->set_label_string(28,
"C3AH6_nucleation_probability");
1279 F->set_label_string(29,
"C3AH6_scale_factor");
1280 F->set_label_string(30,
"FH3_nucleation_probability");
1281 F->set_label_string(31,
"FH3_scale_factor");
1282 F->set_label_string(32,
"Microstructure_size");
1283 F->set_label_string(33,
"Adiabatic_conditions");
1284 F->set_label_string(34,
"Vol_cement_clinker_gypsum");
1285 F->set_label_string(35,
"Vol_cement_SCM");
1286 F->set_label_string(36,
"Vol_water");
1287 F->set_label_string(37,
"Vol_FA");
1288 F->set_label_string(38,
"Vol_CA");
1289 F->set_label_string(39,
"Vol_inert_filler");
1290 F->set_label_string(40,
"Vol_entrained_entrapped_air");
1291 F->set_label_string(41,
"Grain_average_FA");
1292 F->set_label_string(42,
"Grain_average_CA");
1293 F->set_label_string(43,
"ITZ_thickness");
1294 F->set_label_string(44,
"ITZ_Young_red");
1295 F->set_label_string(45,
"Young_SCM");
1296 F->set_label_string(46,
"Poisson_SCM");
1297 F->set_label_string(47,
"Young_FA");
1298 F->set_label_string(48,
"Poisson_FA");
1299 F->set_label_string(49,
"Young_CA");
1300 F->set_label_string(50,
"Poisson_CA");
1301 F->set_label_string(51,
"Young_inert");
1302 F->set_label_string(52,
"Poisson_inert");
1303 F->set_label_string(53,
"Calculate_elastic_homogenization");
1358 F->set_section_string(0,
"CEMHYD_generate_particles");
1360 F->check_requirements();
1362 if ( F->error_in_requirements() ) {
1363 printf(
"Cemhyd input file %s is not complete (file %s, line %d)\n", inp, __FILE__, __LINE__);
1367 F->get_value(0, (
long & )
iseed);
1372 int errorId =
xmlFile->LoadFile(inp);
1373 if ( errorId != XML_SUCCESS ) {
1374 printf(
"\nError reading XML file %s or nonletter symbols used, error id = %d\n", inp, errorId);
1389 F->get_value(32, (
long & )
SYSSIZE);
1392 printf(
"Can not run small microstructure %d (< 10 voxels a side), file %s, line %d)\n",
SYSSIZE, __FILE__, __LINE__);
1423 F->get_value(13, (
long & )read_micr);
1426 if ( !read_micr && generateMicrostructure == 1 ) {
1432 printf(
"MONOPHASE microstructure created\n");
1447double CemhydMatStatus :: ran1(
int *idum)
1456 if ( ( * idum <= 0 ) || (
iy == 0 ) ) {
1457 * idum = ( -* idum > * idum ) ? -* idum : * idum;
1458 for ( j =
NTAB + 7; j >= 0; j-- ) {
1460 * idum =
IA * ( * idum - k *
IQ ) -
IR * k;
1474 * idum =
IA * ( * idum - k *
IQ ) -
IR * k;
1479 j = ( int ) (
iy *
NDIV );
1487void CemhydMatStatus :: addagg()
1496 printf(
"Enter thickness of aggregate to place (an even integer) \n");
1498 F->get_next_line_in_section(0, (
long & )
aggsize);
1515 for ( ix = agglo; ix <= agghi; ix++ ) {
1517 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
1534int CemhydMatStatus :: chksph(
int xin,
int yin,
int zin,
int radd,
int wflg,
int phasein,
int phase2)
1538 int nofits, xp, yp, zp, i, j, k;
1539 float dist, xdist, ydist, zdist, ftmp;
1543 for ( i = xin - radd; ( ( i <= xin + radd ) && ( nofits == 0 ) ); i++ ) {
1552 ftmp = ( float ) ( i - xin );
1553 xdist = ftmp * ftmp;
1554 for ( j = yin - radd; ( ( j <= yin + radd ) && ( nofits == 0 ) ); j++ ) {
1563 ftmp = ( float ) ( j - yin );
1564 ydist = ftmp * ftmp;
1565 for ( k = zin - radd; ( ( k <= zin + radd ) && ( nofits == 0 ) ); k++ ) {
1574 ftmp = ( float ) ( k - zin );
1575 zdist = ftmp * ftmp;
1578 dist = sqrt(xdist + ydist + zdist);
1579 if ( ( dist - 0.5 ) <= (
float ) radd ) {
1582 cement [ xp ] [ yp ] [ zp ] = phasein;
1583 cemreal [ xp ] [ yp ] [ zp ] = phase2;
1586 else if ( ( wflg == 1 ) && (
cement [ xp ] [ yp ] [ zp ] !=
POROSITY ) ) {
1592 if ( ( wflg == 1 ) && ( ( fabs( xp - ( (
float ) (
SYSSIZE + 1 ) / 2.0 ) ) ) < ( (
float )
aggsize / 2.0 ) ) ) {
1611int CemhydMatStatus :: gsphere(
int numgen,
long int *numeach,
int *sizeeach,
int *pheach)
1615 int count, x, y, z, radius, ig, tries, phnow;
1617 float testgyp, typegyp;
1620 for ( ig = 0; ig < numgen; ig++ ) {
1621 phnow = pheach [ ig ];
1622 radius = sizeeach [ ig ];
1624 for ( jg = 1; jg <= numeach [ ig ]; jg++ ) {
1643 printf(
"Could not place sphere %d after %ld random attempts \n",
npart,
MAXTRIES);
1644 printf(
"Skipping this microstructure parameters\n");
1645 for ( i = 1; i <=
npart; i++ ) {
1651 }
while (
count != 0 );
1656 printf(
"Too many spheres being generated \n");
1657 printf(
"User needs to increase value of NPARTC at top of C-code\n");
1658 printf(
"Skipping this microstructure parameters\n");
1713 for ( i = 1; i <=
npart; i++ ) {
1723int CemhydMatStatus :: create()
1728 int *sphrad, *sphase;
1735 sphnum =
new long int [
NUMSIZES ];
1740 printf(
"Enter number of different size spheres to use(max. is %d) \n",
NUMSIZES);
1743 F->get_next_line_in_section(0, (
long & )numsize);
1750 printf(
"%d \n", numsize);
1752 }
while ( ( numsize >
NUMSIZES ) || ( numsize < 0 ) );
1756 printf(
"Enter dispersion factor (separation distance in pixels) for spheres (0-2) \n");
1757 printf(
"0 corresponds to totally random placement \n");
1760 F->get_next_line_in_section(0, (
long & )
dispdist);
1773 printf(
"Enter probability for gypsum particles on a random particle basis (0.0-1.0) \n");
1776 F->get_next_line_in_section(0,
probgyp);
1789 printf(
"Enter probabilities for hemihydrate and anhydrite forms of gypsum (0.0-1.0) \n");
1792 F->get_next_line_in_section(0,
probhem);
1793 F->get_next_line_in_section(0,
probanh);
1805 if ( ( numsize > 0 ) && ( numsize < (
NUMSIZES + 1 ) ) ) {
1807 printf(
"Enter number, radius, and phase ID for each sphere class (largest radius 1st) \n");
1808 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);
1811 for ( isph = 0; isph < numsize; isph++ ) {
1813 printf(
"Enter number of spheres of class %d \n", isph + 1);
1816 F->get_next_line_in_section(0, inval1);
1824 printf(
"%ld \n", inval1);
1826 sphnum [ isph ] = inval1;
1830 printf(
"Enter radius of spheres of class %d \n", isph + 1);
1831 printf(
"(Integer <=%d please) \n",
SYSSIZE / 3);
1834 F->get_next_line_in_section(0, (
long & )inval);
1840 if ( inval > (
SYSSIZE / 3 ) ) {
1841 printf(
"Given radius %d exceeded maximum radius of %d, terminating\n", inval,
SYSSIZE / 3);
1846 printf(
"%d \n", inval);
1850 sphrad [ isph ] = inval;
1853 printf(
"Enter phase of spheres of class %d \n", isph + 1);
1856 F->get_next_line_in_section(0, (
long & )inval);
1863 printf(
"%d \n", inval);
1867 sphase [ isph ] = inval;
1868 if ( inval ==
CEMID ) {
1877 if (
gsphere(numsize, sphnum, sphrad, sphase) == 1 ) {
1896void CemhydMatStatus :: drawfloc(
int xin,
int yin,
int zin,
int radd,
int phasein,
int phase2)
1900 int xp, yp, zp, i, j, k;
1901 float dist, xdist, ydist, zdist, ftmp;
1904 for ( i = xin - radd; ( i <= xin + radd ); i++ ) {
1913 ftmp = ( float ) ( i - xin );
1914 xdist = ftmp * ftmp;
1915 for ( j = yin - radd; ( j <= yin + radd ); j++ ) {
1924 ftmp = ( float ) ( j - yin );
1925 ydist = ftmp * ftmp;
1926 for ( k = zin - radd; ( k <= zin + radd ); k++ ) {
1935 ftmp = ( float ) ( k - zin );
1936 zdist = ftmp * ftmp;
1939 dist = sqrt(xdist + ydist + zdist);
1940 if ( ( dist - 0.5 ) <= (
float ) radd ) {
1942 cement [ xp ] [ yp ] [ zp ] = phasein;
1943 cemreal [ xp ] [ yp ] [ zp ] = phase2;
1954int CemhydMatStatus :: chkfloc(
int xin,
int yin,
int zin,
int radd)
1958 int nofits, xp, yp, zp, i, j, k;
1959 float dist, xdist, ydist, zdist, ftmp;
1964 for ( i = xin - radd; ( ( i <= xin + radd ) && ( nofits == 0 ) ); i++ ) {
1973 ftmp = ( float ) ( i - xin );
1974 xdist = ftmp * ftmp;
1975 for ( j = yin - radd; ( ( j <= yin + radd ) && ( nofits == 0 ) ); j++ ) {
1984 ftmp = ( float ) ( j - yin );
1985 ydist = ftmp * ftmp;
1986 for ( k = zin - radd; ( ( k <= zin + radd ) && ( nofits == 0 ) ); k++ ) {
1995 ftmp = ( float ) ( k - zin );
1996 zdist = ftmp * ftmp;
1999 dist = sqrt(xdist + ydist + zdist);
2000 if ( ( dist - 0.5 ) <= (
float ) radd ) {
2003 nofits =
cement [ xp ] [ yp ] [ zp ];
2008 if ( ( fabs( xp - ( (
float ) (
SYSSIZE + 1 ) / 2.0 ) ) ) < ( (
float )
aggsize / 2.0 ) ) {
2021void CemhydMatStatus :: makefloc()
2025 int partdo, numfloc;
2027 int nleft = 0, ckall;
2028 int xm, ym, zm, moveran;
2029 int xp, yp, zp, rp, clushit, valkeep;
2032 struct cluster *parttmp, *partpoint, *partkeep = NULL;
2034 cluspart =
new int [
NPARTC ];
2038 for ( iclus = 1; iclus <=
npart; iclus++ ) {
2039 cluspart [ iclus ] = iclus;
2044 printf(
"Enter number of flocs desired at end of routine (>0) \n");
2047 F->get_next_line_in_section(0, (
long & )numfloc);
2054 printf(
"%d\n", numfloc);
2056 }
while ( numfloc <= 0 );
2058 while ( nstart > numfloc ) {
2062 for ( iclus = 1; iclus <=
npart; iclus++ ) {
2063 if (
clust [ iclus ] == NULL ) {
2068 moveran = ( int ) ( 6. *
ran1(
seed) );
2069 switch ( moveran ) {
2093 partpoint =
clust [ iclus ];
2094 while ( partpoint != NULL ) {
2105 partpoint =
clust [ iclus ];
2106 while ( ( partpoint != NULL ) && ( ckall == 0 ) ) {
2107 xp = partpoint->
x + xm;
2108 yp = partpoint->
y + ym;
2109 zp = partpoint->
z + zm;
2111 ckall =
chkfloc(xp, yp, zp, rp);
2117 partpoint =
clust [ iclus ];
2118 while ( partpoint != NULL ) {
2119 xp = partpoint->
x + xm;
2120 yp = partpoint->
y + ym;
2121 zp = partpoint->
z + zm;
2124 partdo = partpoint->
partid;
2125 drawfloc(xp, yp, zp, rp, partdo +
CEM - 1, valkeep);
2135 partpoint =
clust [ iclus ];
2137 while ( partpoint != NULL ) {
2143 partdo = partpoint->
partid;
2144 drawfloc(xp, yp, zp, rp, partdo +
CEM - 1, valkeep);
2145 partkeep = partpoint;
2150 if ( ckall !=
AGG ) {
2151 clushit = cluspart [ ckall -
CEM + 1 ];
2153 parttmp =
clust [ clushit ];
2156 while ( parttmp != NULL ) {
2157 cluspart [ parttmp->
partid ] = iclus;
2164 clust [ clushit ] = NULL;
2171 printf(
"Number left was %d but number of clusters is %d \n", nleft, nstart);
2182void CemhydMatStatus :: measure()
2186 long int npor, nc2s, ngyp, ncem, nagg, npozz, ninert, nflyash, nanh, nhem, ncaco3, nslag;
2204 for ( i = 1; i <=
SYSSIZE; i++ ) {
2205 for ( j = 1; j <=
SYSSIZE; j++ ) {
2206 for ( k = 1; k <=
SYSSIZE; k++ ) {
2207 valph =
cemreal [ i ] [ j ] [ k ];
2210 }
else if ( valph ==
CEMID ) {
2212 }
else if ( valph ==
C2SID ) {
2214 }
else if ( valph ==
GYPID ) {
2220 }
else if ( valph ==
AGG ) {
2222 }
else if ( valph ==
POZZID ) {
2224 }
else if ( valph ==
SLAGID ) {
2226 }
else if ( valph ==
INERTID ) {
2228 }
else if ( valph ==
FLYASH ) {
2230 }
else if ( valph ==
CACO3 ) {
2238 printf(
"\n Phase counts are: \n");
2239 printf(
"Porosity= %ld \n", npor);
2240 printf(
"Cement= %ld \n", ncem);
2241 printf(
"C2S= %ld \n", nc2s);
2242 printf(
"Gypsum= %ld \n", ngyp);
2243 printf(
"Anhydrite= %ld \n", nanh);
2244 printf(
"Hemihydrate= %ld \n", nhem);
2245 printf(
"Pozzolan= %ld \n", npozz);
2246 printf(
"Inert= %ld \n", ninert);
2247 printf(
"Slag= %ld \n", nslag);
2248 printf(
"CaCO3= %ld \n", ncaco3);
2249 printf(
"Fly Ash= %ld \n", nflyash);
2250 printf(
"Aggregate= %ld \n", nagg);
2256void CemhydMatStatus :: measagg()
2260 int phase [ 40 ], ptot;
2261 int icnt, ixlo, ixhi,
iy, iz, phid, idist;
2265 aggfile = fopen(
"agglist.out",
"w");
2266 printf(
"Distance Porosity Cement C2S Gypsum Anhydrite Hemihydrate Pozzolan Inert Slag CaCO3 Fly Ash\n");
2267 fprintf(aggfile,
"Distance Porosity Cement C2S Gypsum Anhydrite Hemihydrate Pozzolan Inert Slag CaCO3 Fly Ash\n");
2277 for ( icnt = 0; icnt < 39; icnt++ ) {
2285 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2289 phase [ phid ] += 1;
2295 phase [ phid ] += 1;
2301 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 ]);
2302 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 ]);
2312void CemhydMatStatus :: connect()
2316 long int inew, ntop, nthrough, ncur, nnew, ntot;
2317 int i, j, k, nmatx [ 29000 ], nmaty [ 29000 ], nmatz [ 29000 ];
2318 int xcn, ycn, zcn, npix, x1, y1, z1, igood, nnewx [ 29000 ], nnewy [ 29000 ], nnewz [ 29000 ];
2322 printf(
"Enter phase to analyze 0) pores 1) Cement \n");
2324 F->get_next_line_in_section(0, (
long & )npix);
2330 printf(
"%d \n", npix);
2331 }
while ( ( npix != 0 ) && ( npix != 1 ) );
2341 for ( i = 1; i <=
SYSSIZE; i++ ) {
2342 for ( j = 1; j <=
SYSSIZE; j++ ) {
2346 if ( ( (
cement [ i ] [ j ] [ k ] == npix ) && ( (
cement [ i ] [ j ] [
SYSSIZE ] == npix ) ||
2362 for ( inew = 1; inew <= ncur; inew++ ) {
2363 xcn = nmatx [ inew ];
2364 ycn = nmaty [ inew ];
2365 zcn = nmatz [ inew ];
2368 for ( jnew = 1; jnew <= 6; jnew++ ) {
2377 }
else if ( jnew == 2 ) {
2382 }
else if ( jnew == 3 ) {
2387 }
else if ( jnew == 4 ) {
2392 }
else if ( jnew == 5 ) {
2394 }
else if ( jnew == 6 ) {
2399 if ( ( z1 >= 1 ) && ( z1 <=
SYSSIZE ) ) {
2400 if ( (
cement [ x1 ] [ y1 ] [ z1 ] == npix ) || ( (
cement [ x1 ] [ y1 ] [ z1 ] >=
CEM ) &&
2401 (
cement [ x1 ] [ y1 ] [ z1 ] <
BURNTG ) && ( npix == 1 ) ) ) {
2405 if ( nnew >= 29000 ) {
2406 printf(
"error in size of nnew \n");
2409 nnewx [ nnew ] = x1;
2410 nnewy [ nnew ] = y1;
2411 nnewz [ nnew ] = z1;
2424 for ( icur = 1; icur <= ncur; icur++ ) {
2425 nmatx [ icur ] = nnewx [ icur ];
2426 nmaty [ icur ] = nnewy [ icur ];
2427 nmatz [ icur ] = nnewz [ icur ];
2430 }
while ( nnew > 0 );
2440 printf(
"Phase ID= %d \n", npix);
2441 printf(
"Number accessible from top= %ld \n", ntop);
2442 printf(
"Number contained in through pathways= %ld \n", nthrough);
2445 for ( i = 1; i <=
SYSSIZE; i++ ) {
2446 for ( j = 1; j <=
SYSSIZE; j++ ) {
2447 for ( k = 1; k <=
SYSSIZE; k++ ) {
2458void CemhydMatStatus :: outmic()
2462 FILE *outfile, *partfile;
2463 char filen [ 80 ], filepart [ 80 ];
2464 int ix,
iy, iz, valout;
2467 printf(
"Enter name of file to save microstructure to \n");
2470 F->get_next_line_in_section(0, filen);
2476 printf(
"%s\n", filen);
2478 outfile = fopen(filen,
"w");
2481 printf(
"Enter name of file to save particle IDs to \n");
2484 F->get_next_line_in_section(0, filepart);
2491 printf(
"%s\n", filepart);
2494 partfile = fopen(filepart,
"w");
2496 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2498 for ( ix = 1; ix <=
SYSSIZE; ix++ ) {
2500 fprintf(outfile,
"%1d\n", valout);
2501 valout =
cement [ ix ] [
iy ] [ iz ];
2506 fprintf(partfile,
"%d\n", valout);
2516int CemhydMatStatus :: genpartnew()
2586 for ( ig = 1; ig <=
SYSSIZE; ig++ ) {
2587 for ( jg = 1; jg <=
SYSSIZE; jg++ ) {
2588 for ( kg = 1; kg <=
SYSSIZE; kg++ ) {
2598 printf(
" \n Input User Choice \n");
2599 printf(
"1) Exit \n");
2600 printf(
"2) Add spherical particles (cement,gypsum, pozzolans, etc.) to microstructure \n");
2601 printf(
"3) Flocculate system by reducing number of particle clusters \n");
2602 printf(
"4) Measure global phase fractions \n");
2603 printf(
"5) Add an aggregate to the microstructure \n");
2604 printf(
"6) Measure single phase connectivity (pores or solids) \n");
2605 printf(
"7) Measure phase fractions vs. distance from aggregate surface \n");
2606 printf(
"8) Output current microstructure to file \n");
2609 F->get_next_line_in_section(0, (
long & )userc);
2644 printf(
"No aggregate present. \n");
2654 }
while ( userc != 1 );
2657 for ( ig = 0; ig <
SYSSIZE; ig++ ) {
2658 for ( jg = 0; jg <
SYSSIZE; jg++ ) {
2659 for ( kg = 0; kg <
SYSSIZE; kg++ ) {
2660 micpart [ ig ] [ jg ] [ kg ] =
cement [ ig + 1 ] [ jg + 1 ] [ kg + 1 ];
2679int CemhydMatStatus :: maketemp(
int size)
2681 int icirc, xval, yval, zval;
2687 for ( xval = ( -size ); xval <= size; xval++ ) {
2688 xtmp = ( float ) ( xval * xval );
2689 for ( yval = ( -size ); yval <= size; yval++ ) {
2690 ytmp = ( float ) ( yval * yval );
2691 for ( zval = ( -size ); zval <= size; zval++ ) {
2692 dist = sqrt( xtmp + ytmp + (
float ) ( zval * zval ) );
2693 if ( dist <= ( (
float ) size + 0.5 ) ) {
2696 printf(
"Too many elements in sphere \n");
2697 printf(
"Must change value of MAXSPH parameter \n");
2698 printf(
"Currently set at %ld \n",
MAXSPH);
2702 xsph [ icirc ] = xval;
2703 ysph [ icirc ] = yval;
2704 zsph [ icirc ] = zval;
2717void CemhydMatStatus :: phcount()
2719 long int npore,
nsolid [ 37 ];
2723 for ( ix = 1; ix < 37; ix++ ) {
2728 for ( ix = 1; ix <=
SYSSIZE; ix++ ) {
2730 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2731 if (
mask [ ix ] [
iy ] [ iz ] == 0 ) {
2740 printf(
"Pores are: %ld \n", npore);
2741 printf(
"Solids are: %ld %ld %ld %ld %ld %ld\n",
nsolid [ 1 ],
nsolid [ 2 ],
2749int CemhydMatStatus :: surfpix(
int xin,
int yin,
int zin)
2751 int npix, ix1, iy1, iz1;
2762 if (
mask [ ix1 ] [ yin ] [ zin ] == 0 ) {
2771 if (
mask [ ix1 ] [ yin ] [ zin ] == 0 ) {
2780 if (
mask [ xin ] [ iy1 ] [ zin ] == 0 ) {
2789 if (
mask [ xin ] [ iy1 ] [ zin ] == 0 ) {
2798 if (
mask [ xin ] [ yin ] [ iz1 ] == 0 ) {
2807 if (
mask [ xin ] [ yin ] [ iz1 ] == 0 ) {
2817float CemhydMatStatus :: rhcalc(
int phin)
2820 long int porc, surfc;
2826 for ( ix = 1; ix <=
SYSSIZE; ix++ ) {
2828 for ( iz = 1; iz <=
SYSSIZE; iz++ ) {
2829 if (
mask [ ix ] [
iy ] [ iz ] == phin ) {
2837 printf(
"Phase area count is %ld \n", porc);
2838 printf(
"Phase surface count is %ld \n", surfc);
2839 rhval = ( float ) porc * 6. / ( 4. * (
float ) surfc );
2840 printf(
"Hydraulic radius is %f \n", rhval);
2848int CemhydMatStatus :: countem(
int xp,
int yp,
int zp,
int phin)
2854 for ( ic = 1; ic <=
nsph; ic++ ) {
2855 xc = xp +
xsph [ ic ];
2856 yc = yp +
ysph [ ic ];
2857 zc = zp +
zsph [ ic ];
2877 if ( ( xc != xp ) || ( yc != yp ) || ( zc != zp ) ) {
2878 if ( (
mask [ xc ] [ yc ] [ zc ] == phin ) || (
mask [ xc ] [ yc ] [ zc ] == 0 ) ) {
2891void CemhydMatStatus :: sysinit(
int ph1,
int ph2)
2893 int count, xl, yl, zl;
2897 for ( xl = 1; xl <=
SYSSIZE; xl++ ) {
2898 for ( yl = 1; yl <=
SYSSIZE; yl++ ) {
2899 for ( zl = 1; zl <=
SYSSIZE; zl++ ) {
2903 if (
mask [ xl ] [ yl ] [ zl ] == ph1 ) {
2909 if (
mask [ xl ] [ yl ] [ zl ] == ph2 ) {
2914 printf(
"Error count is %d \n",
count);
2915 printf(
"xl %d yl %d zl %d \n", xl, yl, zl);
2920 if ( (
count >= 0 ) && (
mask [ xl ] [ yl ] [ zl ] == ph1 ) ) {
2927 if ( (
count >= 0 ) && (
mask [ xl ] [ yl ] [ zl ] == ph2 ) ) {
2943void CemhydMatStatus :: sysscan(
int ph1,
int ph2)
2945 int xd, yd, zd, curvval;
2948 for ( xd = 1; xd <=
SYSSIZE; xd++ ) {
2949 for ( yd = 1; yd <=
SYSSIZE; yd++ ) {
2950 for ( zd = 1; zd <=
SYSSIZE; zd++ ) {
2951 curvval =
curvature [ xd ] [ yd ] [ zd ];
2953 if (
mask [ xd ] [ yd ] [ zd ] == ph2 ) {
2954 nair [ curvval ] += 1;
2955 }
else if (
mask [ xd ] [ yd ] [ zd ] == ph1 ) {
2968int CemhydMatStatus :: procsol(
int nsearch)
2970 int valfound, i, stop;
2975 valfound =
nsph - 1;
2978 for ( i = (
nsph - 1 ); ( ( i >= 0 ) && ( stop == 0 ) ); i-- ) {
2980 if ( nsofar > nsearch ) {
2986 return ( valfound );
2995int CemhydMatStatus :: procair(
int nsearch)
2997 int valfound, i, stop;
3005 for ( i = 0; ( ( i <
nsph ) && ( stop == 0 ) ); i++ ) {
3006 nsofar +=
nair [ i ];
3007 if ( nsofar > nsearch ) {
3013 return ( valfound );
3020int CemhydMatStatus :: movepix(
int ntomove,
int ph1,
int ph2)
3022 int count1, count2, ntot, countc, i, xp, yp, zp;
3023 int cmin, cmax, cfg;
3025 long int nsolc, nairc, nsum, nsolm, nairm, nst1, nst2, next1, next2;
3026 float pck, plsol, plair;
3034 for ( i =
nsph; i > count1; i-- ) {
3035 if ( (
nsolid [ i ] > 0 ) && ( cfg == 0 ) ) {
3044 plsol = ( float ) ( ntomove - nsum ) / ( float )
nsolid [ count1 ];
3045 next1 = ntomove - nsum;
3046 nst1 =
nsolid [ count1 ];
3052 for ( i = 0; i < count2; i++ ) {
3053 if ( (
nair [ i ] > 0 ) && ( cfg == 0 ) ) {
3062 plair = ( float ) ( ntomove - nsum ) / ( float )
nair [ count2 ];
3063 next2 = ntomove - nsum;
3064 nst2 =
nair [ count2 ];
3068 if ( cmin >= cmax ) {
3070 printf(
"Stopping - at equilibrium \n");
3071 printf(
"cmin- %d cmax- %d \n", cmin, cmax);
3083 for ( xp = 1; xp <=
SYSSIZE; xp++ ) {
3084 for ( yp = 1; yp <=
SYSSIZE; yp++ ) {
3085 for ( zp = 1; zp <=
SYSSIZE; zp++ ) {
3086 countc =
curvature [ xp ] [ yp ] [ zp ];
3088 if (
mask [ xp ] [ yp ] [ zp ] == ph1 ) {
3089 if ( countc > count1 ) {
3091 mask [ xp ] [ yp ] [ zp ] = ph2;
3095 nair [ countc ] += 1;
3100 if ( countc == count1 ) {
3104 if ( ( pck < 0 ) || ( pck > 1.0 ) ) {
3108 if ( ( ( pck < plsol ) && ( nsolc < next1 ) ) || ( ( nst1 - nsolm ) < ( next1 - nsolc ) ) ) {
3111 mask [ xp ] [ yp ] [ zp ] = ph2;
3115 nair [ count1 ] += 1;
3122 else if (
mask [ xp ] [ yp ] [ zp ] == ph2 ) {
3123 if ( countc < count2 ) {
3125 mask [ xp ] [ yp ] [ zp ] = ph1;
3128 nair [ countc ] -= 1;
3132 if ( countc == count2 ) {
3135 if ( ( pck < 0 ) || ( pck > 1.0 ) ) {
3139 if ( ( ( pck < plair ) && ( nairc < next2 ) ) || ( ( nst2 - nairm ) < ( next2 - nairc ) ) ) {
3142 mask [ xp ] [ yp ] [ zp ] = ph1;
3145 nair [ count2 ] -= 1;
3159 printf(
"ntot is %d \n", ntot);
3166void CemhydMatStatus :: sinter3d(
int ph1id,
int ph2id,
float rhtarget)
3168 int natonce, i, rade, j, rflag;
3170 long int curvsum1, curvsum2, pixsum1, pixsum2;
3171 float rhnow, avecurv1, avecurv2;
3174 for ( i = 0; i <= 499; i++ ) {
3185 printf(
"nsph is %d \n",
nsph);
3194 while ( ( rhnow < rhtarget ) && ( i <
MAXCYC_SEAL ) ) {
3195 printf(
"Now: %f Target: %f \n", rhnow, rhtarget);
3198 printf(
"Cycle: %d \n", i);
3200 keepgo =
movepix(natonce, ph1id, ph2id);
3202 if ( keepgo == 1 ) {
3211 for ( j = 0; j <=
nsph; j++ ) {
3213 curvsum1 += ( j *
nsolid [ j ] );
3214 pixsum2 +=
nair [ j ];
3215 curvsum2 += ( j *
nair [ j ] );
3218 avecurv1 = ( float ) curvsum1 / (
float ) pixsum1;
3219 avecurv2 = ( float ) curvsum2 / (
float ) pixsum2;
3220 printf(
"Ave. solid curvature: %f \n", avecurv1);
3221 printf(
"Ave. air curvature: %f \n", avecurv2);
3226void CemhydMatStatus :: stat3d()
3228 int valin, ix,
iy, iz;
3229 int ix1, iy1, iz1, k;
3230 long int voltot, surftot;
3232 for ( ix = 0; ix <= 42; ix++ ) {
3237 for ( iz = 1; iz <=
SYSIZE; iz++ ) {
3239 for ( ix = 1; ix <=
SYSIZE; ix++ ) {
3240 valin =
mask [ ix ] [
iy ] [ iz ];
3247 for ( iz = 1; iz <=
SYSIZE; iz++ ) {
3249 for ( ix = 1; ix <=
SYSIZE; ix++ ) {
3250 if (
mask [ ix ] [
iy ] [ iz ] != 0 ) {
3251 valin =
mask [ ix ] [
iy ] [ iz ];
3253 for ( k = 1; k <= 6; k++ ) {
3313 if ( ( ix1 < 1 ) || ( iy1 < 1 ) || ( iz1 < 1 ) || ( ix1 >
SYSIZE ) || ( iy1 >
SYSIZE ) || ( iz1 >
SYSIZE ) ) {
3314 printf(
"%d %d %d \n", ix1, iy1, iz1);
3318 if (
mask [ ix1 ] [ iy1 ] [ iz1 ] == 0 ) {
3328 printf(
"Phase Volume Surface Volume Surface \n");
3329 printf(
" ID count count fraction fraction \n");
3338 for ( k = 1; k <= 4; k++ ) {
3339 printf(
" %d %8ld %8ld %.5f %.5f\n", k,
volume [ k ],
surface [ k ],
3340 (
float )
volume [ k ] / (
float ) voltot, (
float )
surface [ k ] / (
float ) surftot);
3343 printf(
"Total %8ld %8ld\n\n\n", voltot, surftot);
3345 for ( k = 5; k <= 11; k++ ) {
3351 for ( k = 24; k <= 27; k++ ) {
3361void CemhydMatStatus :: rand3d(
int phasein,
int phaseout,
float xpt)
3363 float s2, ss, sdiff, xtmp, ytmp;
3366 double ***normm, ***res;
3371 float *s, *xr, *
sum;
3373 double t1, t2, x1, x2, u1, u2, xrad, resmax, resmin;
3374 float xtot, filval, radius, sect, sumtot, vcrit;
3375 int valin, r1, r2, i1, i2, i3, i, j, k, j1, k1;
3376 int ido, iii, jjj, ix,
iy, iz, index;
3378 char tempstr [ 256 ];
3385 s =
new float [ 61 ];
3386 xr =
new float [ 61 ];
3387 sum =
new float [ 502 ];
3394 t1 = 2. *
M_PI * u2;
3395 t2 = sqrt( -2. * log(u1) );
3398 normm [ i1 ] [ i2 ] [ i3 ] = x1;
3409 normm [ i1 ] [ i2 ] [ i3 ] = x2;
3423 F->get_next_line_in_section(0, (
long & )ido);
3431 printf(
"Number of points in correlation file is %d \n", ido);
3434 for ( i = 1; i <= ido; i++ ) {
3436 F->get_next_line_in_section(0, tempstr);
3437 sscanf(tempstr,
"%d %f", & valin, & val2);
3445 s [ i ] = ( float ) val2;
3446 xr [ i ] = ( float ) r [ i ];
3453 for ( i = 0; i < 31; i++ ) {
3455 for ( j = 0; j < 31; j++ ) {
3457 for ( k = 0; k < 31; k++ ) {
3458 xtmp = ( float ) ( iii + jjj + k * k );
3459 radius = sqrt(xtmp);
3460 r1 = ( int ) radius + 1;
3462 if ( s [ r1 ] < 0.0 ) {
3463 printf(
"%d and %d %f and %f with xtmp of %f\n", r1, r2, s [ r1 ], s [ r2 ], xtmp);
3468 xrad = radius + 1 - r1;
3469 filval = s [ r1 ] + ( s [ r2 ] - s [ r1 ] ) * xrad;
3470 filter [ i + 1 ] [ j + 1 ] [ k + 1 ] = ( filval - s2 ) / sdiff;
3479 for ( i = 1; i <=
SYSIZE; i++ ) {
3480 for ( j = 1; j <=
SYSIZE; j++ ) {
3481 for ( k = 1; k <=
SYSIZE; k++ ) {
3482 res [ i ] [ j ] [ k ] = 0.0;
3483 if ( (
float )
mask [ i ] [ j ] [ k ] == phasein ) {
3484 for ( ix = 1; ix <= 31; ix++ ) {
3494 for (
iy = 1;
iy <= 31;
iy++ ) {
3504 for ( iz = 1; iz <= 31; iz++ ) {
3514 res [ i ] [ j ] [ k ] += normm [ i1 ] [ j1 ] [ k1 ] * filter [ ix ] [
iy ] [ iz ];
3519 if ( res [ i ] [ j ] [ k ] > resmax ) {
3520 resmax = res [ i ] [ j ] [ k ];
3523 if ( res [ i ] [ j ] [ k ] < resmin ) {
3524 resmin = res [ i ] [ j ] [ k ];
3535 printf(
"%d out of %d\n", i,
SYSIZE);
3543 sect = ( resmax - resmin ) / 500.;
3544 for ( i = 1; i <= 500; i++ ) {
3549 for ( i = 1; i <=
SYSIZE; i++ ) {
3550 for ( j = 1; j <=
SYSIZE; j++ ) {
3551 for ( k = 1; k <=
SYSIZE; k++ ) {
3552 if ( (
float )
mask [ i ] [ j ] [ k ] == phasein ) {
3554 index = 1 + ( int ) ( ( res [ i ] [ j ] [ k ] - resmin ) / sect );
3555 if ( index > 500 ) {
3559 sum [ index ] += 1.0;
3566 sumtot = vcrit = 0.0;
3568 for ( i = 1; ( ( i <= 500 ) && ( done == 0 ) ); i++ ) {
3569 sumtot +=
sum [ i ] / xtot;
3570 if ( sumtot > xpt ) {
3572 vcrit = resmin + ( resmax - resmin ) * ( ytmp - 0.5 ) / 500.;
3578 printf(
"Critical volume fraction is %f\n", vcrit);
3581 for ( k = 1; k <=
SYSIZE; k++ ) {
3582 for ( j = 1; j <=
SYSIZE; j++ ) {
3583 for ( i = 1; i <=
SYSIZE; i++ ) {
3584 if ( (
float )
mask [ i ] [ j ] [ k ] == phasein ) {
3585 if ( res [ i ] [ j ] [ k ] > vcrit ) {
3586 mask [ i ] [ j ] [ k ] = phaseout;
3604void CemhydMatStatus :: distrib3d()
3606 int i, j, k, alumval, alum2, valin;
3608 double volin, volf [ 5 ], surff [ 5 ], rhtest, rdesire;
3612 FILE *outfile_img = NULL, *outfile_id = NULL;
3622 printf(
"%d\n", *
seed);
3661 for ( i = 1; i <= 4; i++ ) {
3663 F->get_next_line_in_section(0, volin);
3686 printf(
"Volume %f\n", volf [ i ]);
3689 surff [ i ] = volin;
3691 printf(
"Surface %f\n", surff [ i ]);
3696 if ( ( infile = fopen(filen,
"r") ) != NULL ) {
3697 for ( k = 1; k <=
SYSIZE; k++ ) {
3698 for ( j = 1; j <=
SYSIZE; j++ ) {
3699 for ( i = 1; i <=
SYSIZE; i++ ) {
3700 if ( fscanf(infile,
"%d", & valin) != 1 ) {
3704 mask [ i ] [ j ] [ k ] = valin;
3712 for ( k = 1; k <=
SYSIZE; k++ ) {
3713 for ( j = 1; j <=
SYSIZE; j++ ) {
3714 for ( i = 1; i <=
SYSIZE; i++ ) {
3715 mask [ i ] [ j ] [ k ] =
cemreal [ i ] [ j ] [ k ];
3724 volin = volf [ 1 ] + volf [ 2 ];
3725 if ( volin < 1.0 ) {
3726 rand3d(1, alumval, volin);
3730 rdesire = ( surff [ 1 ] + surff [ 2 ] ) * (
float ) (
surface [ 1 ] +
surface [ alumval ] );
3731 if ( rdesire != 0.0 ) {
3732 if ( (
int ) rdesire <
surface [ 1 ] ) {
3733 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 1 ] ) / rdesire;
3736 rdesire = ( surff [ 3 ] + surff [ 4 ] ) * (
float ) (
surface [ 1 ] +
surface [ alumval ] );
3737 if ( rdesire != 0.0 ) {
3738 rhtest = ( 6. / 4. ) * (
float ) (
volume [ alumval ] ) / rdesire;
3746 if ( ( volf [ 1 ] + volf [ 2 ] ) > 0.0 ) {
3747 volin = volf [ 1 ] / ( volf [ 1 ] + volf [ 2 ] );
3748 if ( volin < 1.0 ) {
3753 rdesire = ( surff [ 1 ] / ( surff [ 1 ] + surff [ 2 ] ) ) * (
float ) (
surface [ 1 ] +
surface [ 2 ] );
3754 if ( rdesire != 0.0 ) {
3755 if ( (
int ) rdesire <
surface [ 1 ] ) {
3756 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 1 ] ) / rdesire;
3759 rdesire = ( surff [ 2 ] / ( surff [ 1 ] + surff [ 2 ] ) ) * (
float ) (
surface [ 1 ] +
surface [ 2 ] );
3760 if ( rdesire != 0.0 ) {
3761 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 2 ] ) / rdesire;
3770 if ( alumval == 4 ) {
3771 volin = volf [ 4 ] / ( volf [ 4 ] + volf [ 3 ] );
3774 volin = volf [ 3 ] / ( volf [ 4 ] + volf [ 3 ] );
3778 if ( volin < 1.0 ) {
3779 rand3d(alumval, alum2, volin);
3783 if ( alumval == 4 ) {
3784 rdesire = ( surff [ 4 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * (
float ) (
surface [ 3 ] +
surface [ 4 ] );
3785 if ( rdesire != 0.0 ) {
3786 if ( (
int ) rdesire <
surface [ 4 ] ) {
3787 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 4 ] ) / rdesire;
3790 rdesire = ( surff [ 3 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * (
float ) (
surface [ 3 ] +
surface [ 4 ] );
3791 if ( rdesire != 0.0 ) {
3792 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 3 ] ) / rdesire;
3798 rdesire = ( surff [ 3 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * (
float ) (
surface [ 3 ] +
surface [ 4 ] );
3799 if ( rdesire != 0.0 ) {
3800 if ( (
int ) rdesire <
surface [ 3 ] ) {
3801 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 3 ] ) / rdesire;
3804 rdesire = ( surff [ 4 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * (
float ) (
surface [ 3 ] +
surface [ 4 ] );
3805 if ( rdesire != 0.0 ) {
3806 rhtest = ( 6. / 4. ) * (
float ) (
volume [ 4 ] ) / rdesire;
3818 F->get_value(14, (
long & )output_img);
3827 F->get_value(15, filen);
3832 if ( ( outfile_img = fopen(filen,
"w") ) == NULL ) {
3833 printf(
"Output img file %s can not be created (file %s, line %d)\n", filen, __FILE__, __LINE__);
3838 F->get_value(16, filen);
3843 if ( ( outfile_id = fopen(filen,
"w") ) == NULL ) {
3844 printf(
"Output id file %s can not be created (file %s, line %d)\n", filen, __FILE__, __LINE__);
3849 for ( k = 0; k <
SYSIZE; k++ ) {
3850 for ( j = 0; j <
SYSIZE; j++ ) {
3851 for ( i = 0; i <
SYSIZE; i++ ) {
3852 mic [ i ] [ j ] [ k ] =
mask [ i + 1 ] [ j + 1 ] [ k + 1 ];
3853 micorig [ i ] [ j ] [ k ] =
mic [ i ] [ j ] [ k ];
3855 fprintf(outfile_img,
"%d\n",
mic [ i ] [ j ] [ k ]);
3856 fprintf(outfile_id,
"%ld\n",
micpart [ i ] [ j ] [ k ]);
3863 fclose(outfile_img);
3880void CemhydMatStatus :: init()
3883 double slagin, CHperslag;
3886 for ( i = 0; i <=
EMPTYP; i++ ) {
4220 if ( ( slagfile = fopen(
"slagchar.dat",
"r") ) == NULL ) {
4225 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4229 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4233 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", &
slagcasi) != 1 ) {
4257 if ( fscanf(slagfile,
"%lf", &
slaghydcasi) != 1 ) {
4261 if ( fscanf(slagfile,
"%lf", &
siperslag) != 1 ) {
4265 if ( fscanf(slagfile,
"%lf", & slagin) != 1 ) {
4270 if ( fscanf(slagfile,
"%lf", &
slagc3a) != 1 ) {
4274 if ( fscanf(slagfile,
"%lf", &
slagreact) != 1 ) {
4298 if ( CHperslag < 0.0 ) {
4309 printf(
"Error in range of C3A/slag value... reset to 1.0 \n");
4317int CemhydMatStatus :: chckedge(
int xck,
int yck,
int zck)
4319 int edgeback, x2, y2, z2;
4326 for ( ip = 0; ( ( ip <
NEIGHBORS ) && ( edgeback == 0 ) ); ip++ ) {
4327 x2 = xck +
xoff [ ip ];
4328 y2 = yck +
yoff [ ip ];
4329 z2 = zck +
zoff [ ip ];
4359 return ( edgeback );
4366void CemhydMatStatus :: passone(
int low,
int high,
int cycid,
int cshexflag)
4368 int i, xid, yid, zid, phid, edgef, phread, cshcyc;
4376 for ( i = low; i <= high; i++ ) {
4381 for ( xid = 0; xid <
SYSIZE; xid++ ) {
4382 for ( yid = 0; yid <
SYSIZE; yid++ ) {
4383 for ( zid = 0; zid <
SYSIZE; zid++ ) {
4384 phread =
mic [ xid ] [ yid ] [ zid ];
4386 if ( ( cshexflag == 1 ) && ( phread ==
CSH ) ) {
4387 cshcyc =
cshage [ xid ] [ yid ] [ zid ];
4394 for ( i = low; ( ( i <= high ) && ( phid == 60 ) ); i++ ) {
4395 if (
mic [ xid ] [ yid ] [ zid ] == i ) {
4410 else if ( i ==
C3S ) {
4412 }
else if ( i ==
C2S ) {
4414 }
else if ( i ==
C3A ) {
4416 }
else if ( i ==
C4AF ) {
4418 }
else if ( i ==
GYPSUM ) {
4426 }
else if ( i ==
POZZ ) {
4428 }
else if ( i ==
SLAG ) {
4430 }
else if ( i ==
ETTR ) {
4441 if ( ( cycid != 0 ) && (
soluble [ phid ] == 1 ) ) {
4464int CemhydMatStatus :: loccsh(
int xcur,
int ycur,
int zcur,
int extent)
4466 int effort, tries, xmod, ymod, zmod;
4467 struct ants *antnew;
4472 while ( ( effort == 0 ) && ( tries < 500 ) ) {
4474 xmod = ( -extent ) + (
int ) ( ( 2. * ( float ) extent + 1. ) *
ran1(
seed) );
4475 ymod = ( -extent ) + (
int ) ( ( 2. * ( float ) extent + 1. ) *
ran1(
seed) );
4476 zmod = ( -extent ) + (
int ) ( ( 2. * ( float ) extent + 1. ) *
ran1(
seed) );
4477 if ( xmod > extent ) {
4481 if ( ymod > extent ) {
4485 if ( zmod > extent ) {
4495 }
else if ( xmod < 0 ) {
4501 }
else if ( zmod < 0 ) {
4507 }
else if ( ymod >=
SYSIZE ) {
4511 if (
mic [ xmod ] [ ymod ] [ zmod ] ==
POROSITY ) {
4517 antnew = (
struct ants * ) malloc(
sizeof(
struct ants ) );
4538int CemhydMatStatus :: countbox(
int boxsize,
int qx,
int qy,
int qz)
4540 int nfound, ix,
iy, iz, qxlo, qxhi, qylo, qyhi, qzlo, qzhi;
4541 int hx, hy, hz, boxhalf;
4543 boxhalf = boxsize / 2;
4545 qxlo = qx - boxhalf;
4546 qxhi = qx + boxhalf;
4547 qylo = qy - boxhalf;
4548 qyhi = qy + boxhalf;
4549 qzlo = qz - boxhalf;
4550 qzhi = qz + boxhalf;
4553 for ( ix = qxlo; ix <= qxhi; ix++ ) {
4557 }
else if ( hx >=
SYSIZE ) {
4561 for (
iy = qylo;
iy <= qyhi;
iy++ ) {
4565 }
else if ( hy >=
SYSIZE ) {
4569 for ( iz = qzlo; iz <= qzhi; iz++ ) {
4573 }
else if ( hz >=
SYSIZE ) {
4578 if ( (
mic [ hx ] [ hy ] [ hz ] <
C3S ) || (
mic [ hx ] [ hy ] [ hz ] >
ABSGYP ) ) {
4594int CemhydMatStatus :: countboxc(
int boxsize,
int qx,
int qy,
int qz)
4596 int nfound, ix,
iy, iz, qxlo, qxhi, qylo, qyhi, qzlo, qzhi;
4597 int hx, hy, hz, boxhalf;
4599 boxhalf = boxsize / 2;
4601 qxlo = qx - boxhalf;
4602 qxhi = qx + boxhalf;
4603 qylo = qy - boxhalf;
4604 qyhi = qy + boxhalf;
4605 qzlo = qz - boxhalf;
4606 qzhi = qz + boxhalf;
4609 for ( ix = qxlo; ix <= qxhi; ix++ ) {
4613 }
else if ( hx >=
SYSIZE ) {
4617 for (
iy = qylo;
iy <= qyhi;
iy++ ) {
4621 }
else if ( hy >=
SYSIZE ) {
4625 for ( iz = qzlo; iz <= qzhi; iz++ ) {
4629 }
else if ( hz >=
SYSIZE ) {
4634 if ( (
mic [ hx ] [ hy ] [ hz ] <
C3S ) || (
mic [ hx ] [ hy ] [ hz ] >
POZZ ) ) {
4652void CemhydMatStatus :: makeinert(
long int ndesire) {
4653 struct togo *headtogo, *tailtogo, *newtogo =
nullptr, *lasttogo, *onetogo;
4655 int px, py, pz, placed, cntpore, cntmax;
4658 headtogo = (
struct togo * ) malloc(
sizeof(
struct togo ) );
4659 headtogo->
x = headtogo->
y = headtogo->
z = ( -1 );
4660 headtogo->
npore = 0;
4663 tailtogo = headtogo;
4667 printf(
"In makeinert with %ld needed elements \n", ndesire);
4672 for ( idesire = 2; idesire <= ndesire; idesire++ ) {
4673 newtogo = (
struct togo * ) malloc(
sizeof(
struct togo ) );
4675 newtogo->
x = newtogo->
y = newtogo->
z = ( -1 );
4682 for ( px = 0; px <
SYSIZE; px++ ) {
4683 for ( py = 0; py <
SYSIZE; py++ ) {
4684 for ( pz = 0; pz <
SYSIZE; pz++ ) {
4687 if ( cntpore > cntmax ) {
4693 if ( cntpore > ( tailtogo->
npore ) ) {
4695 lasttogo = tailtogo;
4696 while ( placed == 0 ) {
4698 if ( newtogo == NULL ) {
4701 if ( cntpore <= ( newtogo->
npore ) ) {
4706 if ( placed == 0 ) {
4711 onetogo = (
struct togo * ) malloc(
sizeof(
struct togo ) );
4715 onetogo->npore = cntpore;
4717 if ( placed == 2 ) {
4718 onetogo->prevtogo = NULL;
4719 onetogo->nexttogo = headtogo;
4724 if ( placed == 1 ) {
4725 onetogo->nexttogo = lasttogo;
4726 onetogo->prevtogo = newtogo;
4732 lasttogo = tailtogo;
4746 for ( idesire = 1; idesire <= ndesire; idesire++ ) {
4750 if ( px != ( -1 ) ) {
4756 lasttogo = headtogo;
4776void CemhydMatStatus :: extslagcsh(
int xpres,
int ypres,
int zpres) {
4777 int check, sump, xchr, ychr, zchr, fchr, i1, action, numnear;
4779 int mstest = 0, mstest2 = 0;
4788 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
4794 sump *=
moveone(& xchr, & ychr, & zchr, & action, sump);
4795 if ( action == 0 ) {
4796 printf(
"Error in value of action in extpozz \n");
4799 check =
mic [ xchr ] [ ychr ] [ zchr ];
4802 if ( xchr != xpres ) {
4807 if ( ychr != ypres ) {
4812 if ( zchr != zpres ) {
4819 if ( (
faces [ xpres ] [ ypres ] [ zpres ] == 0 ) || ( mstest ==
faces [ xpres ] [ ypres ] [ zpres ] ) || ( mstest2 ==
faces [ xpres ] [ ypres ] [ zpres ] ) ) {
4821 faces [ xchr ] [ ychr ] [ zchr ] =
faces [ xpres ] [ ypres ] [ zpres ];
4832 while ( fchr == 0 ) {
4850 check =
mic [ xchr ] [ ychr ] [ zchr ];
4856 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
4869void CemhydMatStatus :: dissolve(
int cycle)
4871 int nc3aext, ncshext, nchext, ngypext, nanhext, plok;
4872 int nsum5, nsum4, nsum3, nsum2, nhemext, nsum6, nc4aext;
4873 int phid, phnew, plnew, cread;
4874 int i, xloop, yloop, zloop, xc, yc;
4878 int placed, cshrand, maxsulfate, msface;
4879 long int ncshgo, nsurf, suminit;
4880 long int xext, nhgd, npchext, nslagc3a = 0;
4881 float pdis, plfh3, fchext, fc3aext, fanhext, mass_now, mass_fa_now, tot_mass, heatfill;
4882 float dfact, dfact1, molesdh2o, h2oinit, heat4, fhemext, fc4aext;
4883 float pconvert, pc3scsh, pc2scsh, calcx, calcy, calcz, tdisfact;
4884 float frafm, frettr, frhyg, frtot, mc3ar, mc4ar, p3init;
4885 struct ants *antadd;
4889 npchext = ncshgo = cshrand = 0;
4895 for ( i = 0; i <=
EMPTYP; i++ ) {
4906 printf(
"Returned from passone \n");
5001 printf(
"Calculated w/c is %.4f\n",
w_to_c);
5002 printf(
"Calculated s/c is %.4f \n",
s_to_c);
5003 printf(
"Calculated heat conversion factor is %f \n",
heat_cf);
5004 printf(
"Calculated mass fractions of water and filler are %.4f and %.4f \n",
5021 printf(
"ctest is %ld\n", ctest);
5028 for ( i = 0; i <=
EMPTYP; i++ ) {
5052 }
else if ( i ==
DIFFCH ) {
5068 }
else if ( i ==
DIFFAS ) {
5086 }
else if ( i ==
C3S ) {
5090 }
else if ( i ==
C2S ) {
5094 }
else if ( i ==
C3A ) {
5099 if ( ( mc3ar + mc4ar ) > 0.0 ) {
5107 frtot = frafm + frettr + frhyg;
5108 if ( frtot > 0.0 ) {
5116 }
else if ( i ==
C4AF ) {
5121 if ( ( mc3ar + mc4ar ) > 0.0 ) {
5128 frtot = frettr + frhyg;
5129 if ( frtot > 0.0 ) {
5161 if ( suminit != 0 ) {
5205 fprintf(
heatfile,
"Cycle time(h) alpha_vol alpha_mass heat4(kJ/kg_solid) Gsratio2 G-s_ratio\n");
5220 printf(
"All water consumed at cycle %d \n",
cyccnt);
5243 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");
5244 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");
5254 for ( i = 0; i <=
HDCSH; i++ ) {
5305 1.42 * (
float )
anhinit + 1.4 * (
float )
heminit + ( (
float )
netbar / 3.30 ) ) ) < 0.25 ) ) ) {
5308 printf(
"Ettringite is soluble beginning at cycle %d \n", cycle);
5359 printf(
"CH dissolution probability goes from %f ",
disprob [
CH ]);
5379 1.42 * (
float )
anhinit + 1.4 * (
float )
heminit ) * 0.05 ) ) ||
5397 if ( maxsulfate > 0 ) {
5722 for ( xloop = 0; xloop <
SYSIZE; xloop++ ) {
5723 for ( yloop = 0; yloop <
SYSIZE; yloop++ ) {
5724 for ( zloop = 0; zloop <
SYSIZE; zloop++ ) {
5725 if (
mic [ xloop ] [ yloop ] [ zloop ] >
OFFSET ) {
5726 phid =
mic [ xloop ] [ yloop ] [ zloop ] -
OFFSET;
5729 if ( ( plnew < 0 ) || ( plnew >=
NEIGHBORS ) ) {
5733 xc = xloop +
xoff [ plnew ];
5734 yc = yloop +
yoff [ plnew ];
5735 zc = zloop +
zoff [ plnew ];
5768 count [ phid ] -= 1;
5770 if ( phid ==
C3AH6 ) {
5775 if ( phid ==
C4AF ) {
5777 if ( ( plfh3 < 0.0 ) || ( plfh3 > 1.0 ) ) {
5783 if ( plfh3 <= 0.5453 ) {
5796 count [ phnew ] += 1;
5797 mic [ xc ] [ yc ] [ zc ] = phnew;
5798 antadd = (
struct ants * ) malloc(
sizeof(
struct ants ) );
5812 if ( ( phid ==
C3S ) || ( phid ==
C2S ) ) {
5814 if ( ( ( phid ==
C2S ) && ( plfh3 <= pc2scsh ) ) || ( plfh3 <= pc3scsh ) ) {
5826 if ( placed != 0 ) {
5835 if ( ( phid ==
C2S ) && ( pc2scsh > 1.0 ) ) {
5837 if ( plfh3 <= ( pc2scsh - 1.0 ) ) {
5849 if ( placed != 0 ) {
5858 mic [ xloop ] [ yloop ] [ zloop ] -=
OFFSET;
5868 if (
mic [ xloop ] [ yloop ] [ zloop ] ==
CSH ) {
5869 if ( (
countbox(3, xloop, yloop, zloop) ) >= 1 ) {
5879 cycnew =
cshage [ xloop ] [ yloop ] [ zloop ];
5881 if ( calcy > 1.0 ) {
5882 calcz = calcy - 1.0;
5884 printf(
"Problem of not creating enough pozzolanic CSH during CSH conversion \n");
5885 printf(
"Current temperature is %f C\n",
temp_cur);
5888 if ( plfh3 <= calcy ) {
5892 mic [ xloop ] [ yloop ] [ zloop ] =
DIFFCH;
5897 antadd = (
struct ants * ) malloc(
sizeof(
struct ants ) );
5922 calcx = ( 19.86 /
molarvcsh [ cycnew ] ) - ( 1. - calcy );
5924 if ( plfh3 < calcx ) {
5933 if (
mic [ xloop ] [ yloop ] [ zloop ] ==
SLAG ) {
5934 if ( (
countbox(3, xloop, yloop, zloop) ) >= 1 ) {
5951 msface = ( int ) ( 3. *
ran1(
seed) + 1. );
5956 faces [ xloop ] [ yloop ] [ zloop ] = msface;
5962 mic [ xloop ] [ yloop ] [ zloop ] =
EMPTYP;
5972 while ( p3init > 1.0 ) {
5978 if ( plfh3 < p3init ) {
5994 if ( ncshgo != 0 ) {
5995 printf(
"CSH dissolved is %ld \n", ncshgo);
5998 if ( npchext > 0 ) {
5999 printf(
"npchext is %ld at cycle %d \n", npchext, cycle);
6006 if ( cshrand != 0 ) {
6007 printf(
"cshrand is %d \n", cshrand);
6021 nchext = ( int ) fchext;
6022 if ( fchext > (
float ) nchext ) {
6024 if ( ( fchext - (
float ) nchext ) > pdis ) {
6043 nc3aext = ( int ) ( fc3aext + nslagc3a );
6044 if ( fc3aext > (
float ) nc3aext ) {
6046 if ( ( fc3aext - (
float ) nc3aext ) > pdis ) {
6052 nc4aext = ( int ) fc4aext;
6053 if ( fc4aext > (
float ) nc4aext ) {
6055 if ( ( fc4aext - (
float ) nc4aext ) > pdis ) {
6068 nanhext = ( int ) fanhext;
6069 if ( fanhext > (
float ) nanhext ) {
6071 if ( ( fanhext - (
float ) nanhext ) > pdis ) {
6083 nhemext = ( int ) fhemext;
6084 if ( fhemext > (
float ) nhemext ) {
6086 if ( ( fhemext - (
float ) nhemext ) > pdis ) {
6099 nsum2 = nchext + ncshext;
6100 nsum3 = nsum2 + nc3aext;
6101 nsum4 = nsum3 + nc4aext;
6102 nsum5 = nsum4 + ngypext;
6103 nsum6 = nsum5 + nhemext;
6107 for ( xext = 1; xext <= ( nsum6 + nanhext ); xext++ ) {
6132 if ( xext > nsum6 ) {
6134 }
else if ( xext > nsum5 ) {
6136 }
else if ( xext > nsum4 ) {
6138 }
else if ( xext > nsum3 ) {
6140 }
else if ( xext > nsum2 ) {
6142 }
else if ( xext > nchext ) {
6146 mic [ xc ] [ yc ] [ zc ] = phid;
6149 antadd = (
struct ants * ) malloc(
sizeof(
struct ants ) );
6165 printf(
"Could not place extra diffusing species (too few POROSITY voxels), continuing (line %d),\n", __LINE__);
6169 }
while ( plok == 0 );
6175 printf(
"Dissolved- %ld %ld %ld %ld %ld %ld %ld %ld %ld %ld %ld %ld\n",
count [
DIFFCSH ],
6198 printf(
"C3AH6 dissolved- %ld with prob. of %f \n", nhgd,
disprob [
C3AH6 ]);
6211void CemhydMatStatus :: addrand(
int randid,
long int nneed)
6215 int success, cpores;
6218 for ( ic = 1; ic <= nneed; ic++ ) {
6220 while ( success == 0 ) {
6237 if ( ( randid !=
CACO3 ) && ( randid !=
INERT ) ) {
6238 mic [ ix ] [
iy ] [ iz ] = randid;
6243 if ( cpores >= 26 ) {
6244 mic [ ix ] [
iy ] [ iz ] = randid;
6255void CemhydMatStatus :: measuresurf()
6257 int sx, sy, sz, jx = 0, jy = 0, jz = 0, faceid;
6259 for ( sx = 0; sx <
SYSIZE; sx++ ) {
6260 for ( sy = 0; sy <
SYSIZE; sy++ ) {
6261 for ( sz = 0; sz <
SYSIZE; sz++ ) {
6263 for ( faceid = 0; faceid < 6; faceid++ ) {
6264 if ( faceid == 1 ) {
6272 }
else if ( faceid == 0 ) {
6274 if ( jx > (
SYSIZE - 1 ) ) {
6280 }
else if ( faceid == 2 ) {
6282 if ( jy > (
SYSIZE - 1 ) ) {
6288 }
else if ( faceid == 3 ) {
6296 }
else if ( faceid == 4 ) {
6298 if ( jz > (
SYSIZE - 1 ) ) {
6304 }
else if ( faceid == 5 ) {
6315 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 ) ) {
6317 if ( (
mic [ jx ] [ jy ] [ jz ] ==
C3S ) || (
mic [ jx ] [ jy ] [ jz ] ==
C2S ) || (
mic [ jx ] [ jy ] [ jz ] ==
C3A ) || (
mic [ jx ] [ jy ] [ jz ] ==
C4AF ) ) {
6328 printf(
"Cement surface count is %ld \n",
scntcement);
6329 printf(
"Total surface count is %ld \n",
scnttotal);
6335 printf(
"Surface fraction is %f \n",
surffract);
6343void CemhydMatStatus :: resaturate()
6346 long int nresat = 0;
6348 for ( sx = 0; sx <
SYSIZE; sx++ ) {
6349 for ( sy = 0; sy <
SYSIZE; sy++ ) {
6350 for ( sz = 0; sz <
SYSIZE; sz++ ) {
6351 if (
mic [ sx ] [ sy ] [ sz ] ==
EMPTYP ) {
6363 printf(
"Number resaturated is %ld \n", nresat);
6369void CemhydMatStatus :: outputImageFileUnperc(
char ***m)
6372 char extension [ 10 ];
6375 prefix = (
char * ) malloc(80);
6383 sprintf(extension,
"%04d",
icyc);
6384 strcpy(prefix,
"unperc/out5.");
6385 strcat(prefix, extension);
6386 strcat(prefix,
".img");
6388 printf(
"Name of output file is %s\n", prefix);
6391 if ( ( imgout = fopen(prefix,
"w") ) == NULL ) {
6392 printf(
"File %s can not be opened\n", prefix);
6397 for (
int dz = 0; dz <
SYSIZE; dz++ ) {
6398 for (
int dy = 0; dy <
SYSIZE; dy++ ) {
6399 for (
int dx = 0; dx <
SYSIZE; dx++ ) {
6400 fprintf(imgout,
"%d\n", m [ dx ] [ dy ] [ dz ]);
6407 printf(
"Unpercolated file %s wrote\n", prefix);
6413void CemhydMatStatus :: readhydrparam()
6415 int fidc3s, fidc2s, fidc3a, fidc4af, fidgyp, fidagg, ffac3a;
6416 int fidhem, fidanh, fidcaco3;
6419 int ix,
iy, iz, phtodo;
6431 printf(
"Seed %d\n", *
seed);
6436 F->get_value(13, (
long & )read_micr);
6445 printf(
"Enter name of file to read initial microstructure from \n");
6448 F->get_value(1, filei);
6456 printf(
"%s\n", filei);
6465 printf(
"Enter IDs in file for C3S, C2S, C3A, C4AF, Gypsum, Hemihydrate, Anhydrite, Aggregate, CaCO3\n");
6480 printf(
"%d %d %d %d %d %d %d %d %d\n", fidc3s, fidc2s, fidc3a, fidc4af, fidgyp, fidhem, fidanh, fidagg, fidcaco3);
6481 printf(
"Enter ID in file for C3A in fly ash (default=35)\n");
6486 printf(
"%d\n", ffac3a);
6491 if ( ( infile = fopen(filei,
"r") ) == NULL ) {
6492 printf(
"CEMHYD3D microstructure file %s not found\n", filei);
6496 for ( iz = 0; iz <
SYSIZE; iz++ ) {
6498 for ( ix = 0; ix <
SYSIZE; ix++ ) {
6500 faces [ ix ] [
iy ] [ iz ] = 0;
6502 if ( fscanf(infile,
"%ld", & valin) == EOF ) {
6503 printf(
"End of file %s reached, terminating\n", filei);
6508 printf(
"Error in the reading at x=%d y=%d z=%d, value %ld\n", ix,
iy, iz, valin);
6512 mic [ ix ] [
iy ] [ iz ] = valin;
6513 if ( valin == fidc3s ) {
6515 }
else if ( valin == fidc2s ) {
6517 }
else if ( ( valin == fidc3a ) || ( valin == ffac3a ) ) {
6519 }
else if ( valin == fidc4af ) {
6521 }
else if ( valin == fidgyp ) {
6523 }
else if ( valin == fidanh ) {
6525 }
else if ( valin == fidhem ) {
6527 }
else if ( valin == fidcaco3 ) {
6529 }
else if ( valin == fidagg ) {
6541 printf(
"Generated microstructure from memory will be used\n");
6550 printf(
"Enter name of file to read particle IDs from \n");
6553 F->get_value(2, filei);
6561 printf(
"%s\n", filei);
6564 if ( ( infile = fopen(filei,
"r") ) == NULL ) {
6565 printf(
"CEMHYD3D microstructure ID file %s not found\n", filei);
6569 for ( iz = 0; iz <
SYSIZE; iz++ ) {
6571 for ( ix = 0; ix <
SYSIZE; ix++ ) {
6572 if ( fscanf(infile,
"%ld", & valin) == EOF ) {
6573 printf(
"End of file %s reached, terminating\n", filei);
6585 printf(
"Generated microstructure ID from memory will be used\n");
6596 printf(
"Enter number of one pixel particles to add (0 to quit) \n");
6602 printf(
"%ld\n", nadd);
6605 while ( nadd > 0 ) {
6607 printf(
"Enter phase to add \n");
6608 printf(
" C3S 1 \n");
6609 printf(
" C2S 2 \n");
6610 printf(
" C3A 3 \n");
6611 printf(
" C4AF 4 \n");
6612 printf(
" GYPSUM 5 \n");
6613 printf(
" HEMIHYD 6 \n");
6614 printf(
" ANHYDRITE 7 \n");
6615 printf(
" POZZ 8 \n");
6616 printf(
" INERT 9 \n");
6617 printf(
" SLAG 10 \n");
6618 printf(
" ASG 11 \n");
6619 printf(
" CAS2 12 \n");
6620 printf(
" CH 13 \n");
6621 printf(
" CSH 14 \n");
6622 printf(
" C3AH6 15 \n");
6623 printf(
" Ettringite 16 \n");
6624 printf(
" Stable Ettringite from C4AF 17 \n");
6625 printf(
" AFM 18 \n");
6626 printf(
" FH3 19 \n");
6627 printf(
" POZZCSH 20 \n");
6628 printf(
" SLAGCSH 21 \n");
6629 printf(
" CACL2 22 \n");
6630 printf(
" Friedels salt 23 \n");
6631 printf(
" Stratlingite 24 \n");
6632 printf(
" Calcium carbonate 26 \n");
6637 printf(
"%d \n", phtodo);
6639 if ( ( phtodo < 0 ) || ( phtodo >
CACO3 ) ) {
6640 printf(
"Error in phase input for one pixel particles \n");
6647 printf(
"Enter number of one pixel particles to add (0 to quit) \n");
6653 printf(
"%ld\n", nadd);
6661 F->get_value(3, (
long & )
sealed);
6671 F->get_value(23, (
long & )
ntimes);
6681 F->get_value(24,
pnucch);
6707 F->get_value(28,
pnuchg);
6733 F->get_value(17, (
long & )
burnfreq);
6746 printf(
"Enter cycle frequency for checking percolation of solids (set) \n");
6752 F->get_value(18, (
long & )
setfreq);
6759 printf(
"Enter cycle frequency for checking hydration of particles \n");
6783 printf(
"Enter apparent activation energy for hydration in kJ/mole \n");
6786 F->get_value(5,
E_act);
6793 printf(
"%f \n",
E_act);
6794 printf(
"Enter apparent activation energy for pozzolanic reactions in kJ/mole \n");
6805 printf(
"Enter apparent activation energy for slag reactions in kJ/mole \n");
6816 printf(
"Enter kinetic factor to convert cycles to time for 25 C \n");
6819 F->get_value(8,
beta);
6826 printf(
"%f \n",
beta);
6827 printf(
"Enter mass fraction of aggregate in concrete \n");
6838 printf(
"Enter kg of cement (+admixtures) in 1 m3\n");
6849 printf(
"Enter heat capacity of aggregate in concrete J/g/C\n");
6852 F->get_value(11,
Cp_agg);
6862 printf(
"Enter heat capacity of cement J/g/C\n");
6875 printf(
"CSH to pozzolanic CSH 0) prohibited or 1) allowed \n");
6881 printf(
"CH precipitation on aggregate surfaces 0) prohibited or 1) allowed \n");
6887 printf(
"Enter number of cycles before executing total resaturation \n");
6893 printf(
"Enter choice for C-S-H geometry 0) random or 1) plates \n");
6899 printf(
"Does pH influence hydration kinetics 0) no or 1) yes \n");
6910 F->get_value(37,
Vol_FA);
6911 F->get_value(38,
Vol_CA);
6966void CemhydMatStatus :: disrealnew_init()
7055 printf(
"After init routine \n");
7062 fileperc = fopen(
"percpore.out",
"w");
7064 fprintf(
fileperc,
"#cycle alpha #through #tot_phase\n");
7073 strcpy(
pHname,
"pHfile.out");
7075 fprintf(
pHfile,
"Cycle time(h) alpha_mass pH sigma [Na+] [K+] [Ca++] [SO4--] activityCa activityOH activitySO4 activityK molesSyngenite\n");
7076 strcpy(fileo,
"out5.img");
7079 strcpy(
ppsname,
"percpore.out");
7080 strcpy(
phrname,
"parthydr.out");
7083 fprintf(
adiafile,
"Cyc Time(h) Temperature Alpha Krate Cp_now Mass_cem kpozz/khyd kslag/khyd DoR_Blend\n");
7085 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");
7086 CSHfile = fopen(
"CSH.out",
"w");
7087 fprintf(
CSHfile,
"cyc alpha [h] HD CSH-total HD/CSH\n");
7088 system(
"mkdir perc 2> /dev/null");
7089 system(
"mkdir unperc 2> /dev/null");
7093 infoperc = fopen(
"perc/info.dat",
"w");
7094 fprintf(
infoperc,
"#Cyc DoH Time[h]\n");
7102 krate = exp( -( 1000. *
E_act / 8.314 ) * ( ( 1. / (
temp_cur + 273.15 ) ) - ( 1. / 298.15 ) ) );
7124void CemhydMatStatus :: disrealnew(
double GiveTemp,
double hydrationTime,
int flag)
7146 while ( flag == 0 ? (
time_cur < hydrationTime ) : counter < flag ) {
7212 krate = exp( -( 1000. *
E_act / 8.314 ) * ( ( 1. / (
temp_cur + 273.15 ) ) - ( 1. / 298.15 ) ) );
7252 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);
7253 fprintf(
heatfile,
"%d %f %f %f %f %f %f \n",
7279 printf(
"Switching to self-desiccating at cycle %d \n",
cyccnt);
7319 double E_paste_SCM_filler_air, nu_paste_SCM_filler_air, E_mortar, nu_mortar;
7326 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 ]);
7367double CemhydMatStatus :: GivePower(
double GiveTemp,
double TargTime)
7370 double castingTime = this->
gp->giveMaterial()->giveCastingTime();
7372 double castingTime = 0.;
7376 if ( TargTime - castingTime <= 0 ) {
7385 printf(
"Cannot go backwards in hydration, TargTime = %f s < LastCallTime = %f s\n", TargTime,
LastCallTime);
7391 if ( !cemhydmat->
nowarnings.
at(4) && ( GiveTemp > 200. ) ) {
7392 printf(
"Temperature exceeds 200 C (file %s, line %d),\n", __FILE__, __LINE__);
7396 if ( GiveTemp > 200. ) {
7397 printf(
"Temperature exceeds 200 C (file %s, line %d),\n", __FILE__, __LINE__);
7406 disrealnew(GiveTemp, ( TargTime - castingTime ) / 3600., 0);
7410 printf(
"time_cur %f [h]\n",
time_cur);
7446 printf(
"PartHeat %f [W/m3 of concrete]\n",
PartHeat);
7455double CemhydMatStatus :: MoveCycles(
double GiveTemp,
int cycles)
7458 double castingTime = 0.;
7460 castingTime = this->
gp->giveMaterial()->giveCastingTime();
7493int CemhydMatStatus :: MoveToDoH(
double GiveTemp,
double DesiredDoH,
int maxcyc)
7498 if ( cycle > maxcyc ) {
7500 printf(
"Target DoH %f (now is %f) was not reached after %d cycles (%s, line %d)\n", DesiredDoH,
alpha_cur, maxcyc, __FILE__, __LINE__);
7515int CemhydMatStatus :: MoveToTime(
double GiveTemp,
double TargTime)
7523double CemhydMatStatus :: GiveTotCemHeat()
7529double CemhydMatStatus :: GiveTotHeat()
7535double CemhydMatStatus :: GiveCp()
7542double CemhydMatStatus :: computeConcreteCapacityBentz()
7544 double capacityPaste, capacityConcrete;
7546 capacityPaste *= ( 1. - 0.26 * ( 1 - pow(2.71828, -2.9 *
alpha_cur) ) );
7549 return capacityConcrete;
7552double CemhydMatStatus :: GiveDensity()
7558double CemhydMatStatus :: GiveDoHLastCyc()
7564double CemhydMatStatus :: GiveDoHActual()
7566 double castingTime = 0;
7573 castingTime = this->
gp->giveMaterial()->giveCastingTime();
7581int CemhydMatStatus :: GiveCycNum()
7587double CemhydMatStatus :: GiveCycTime()
7599int CemhydMatStatus :: burn3d(
int npix,
int d1,
int d2,
int d3)
7603 long int ntop, nthrough, ncur, nnew, ntot, nphc;
7605 int xl, xh, j1, k1, px, py, pz, qx, qy, qz, xcn, ycn, zcn;
7606 int x1, y1, z1, igood;
7607 int *nmatx, *nmaty, *nmatz, *nnewx, *nnewy, *nnewz;
7610 float mass_burn = 0.0, alpha_burn = 0.0;
7612 nmatx =
new int [
SIZE2D ];
7613 nmaty =
new int [
SIZE2D ];
7614 nmatz =
new int [
SIZE2D ];
7615 nnewx =
new int [
SIZE2D ];
7616 nnewy =
new int [
SIZE2D ];
7617 nnewz =
new int [
SIZE2D ];
7633 for ( k = 0; k <
SYSIZE; k++ ) {
7634 for ( j = 0; j <
SYSIZE; j++ ) {
7639 px =
cx(i, j, k, d1, d2, d3);
7640 py =
cy(i, j, k, d1, d2, d3);
7641 pz =
cz(i, j, k, d1, d2, d3);
7642 if (
mic [ px ] [ py ] [ pz ] == npix ) {
7655 for ( inew = 1; inew <= ncur; inew++ ) {
7656 xcn = nmatx [ inew ];
7657 ycn = nmaty [ inew ];
7658 zcn = nmatz [ inew ];
7661 for ( jnew = 1; jnew <= 6; jnew++ ) {
7692 }
else if ( y1 < 0 ) {
7699 }
else if ( z1 < 0 ) {
7704 if ( ( x1 >= 0 ) && ( x1 <
SYSIZE ) ) {
7706 px =
cx(x1, y1, z1, d1, d2, d3);
7707 py =
cy(x1, y1, z1, d1, d2, d3);
7708 pz =
cz(x1, y1, z1, d1, d2, d3);
7709 if (
mic [ px ] [ py ] [ pz ] == npix ) {
7714 printf(
"error in size of nnew \n");
7717 nnewx [ nnew ] = x1;
7718 nnewy [ nnew ] = y1;
7719 nnewz [ nnew ] = z1;
7728 for ( icur = 1; icur <= ncur; icur++ ) {
7729 nmatx [ icur ] = nnewx [ icur ];
7730 nmaty [ icur ] = nnewy [ icur ];
7731 nmatz [ icur ] = nnewz [ icur ];
7734 }
while ( nnew > 0 );
7740 for ( j1 = 0; j1 <
SYSIZE; j1++ ) {
7741 for ( k1 = 0; k1 <
SYSIZE; k1++ ) {
7742 px =
cx(xl, j1, k1, d1, d2, d3);
7743 py =
cy(xl, j1, k1, d1, d2, d3);
7744 pz =
cz(xl, j1, k1, d1, d2, d3);
7745 qx =
cx(xh, j1, k1, d1, d2, d3);
7746 qy =
cy(xh, j1, k1, d1, d2, d3);
7747 qz =
cz(xh, j1, k1, d1, d2, d3);
7748 if ( (
mic [ px ] [ py ] [ pz ] ==
BURNT ) && (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) ) {
7752 if (
mic [ px ] [ py ] [ pz ] ==
BURNT ) {
7753 mic [ px ] [ py ] [ pz ] =
BURNT + 1;
7756 if (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) {
7757 mic [ qx ] [ qy ] [ qz ] =
BURNT + 1;
7770 for ( i = 0; i <
SYSIZE; i++ ) {
7771 for ( j = 0; j <
SYSIZE; j++ ) {
7772 for ( k = 0; k <
SYSIZE; k++ ) {
7773 if (
mic [ i ] [ j ] [ k ] >=
BURNT ) {
7775 mic [ i ] [ j ] [ k ] = npix;
7776 }
else if (
mic [ i ] [ j ] [ k ] == npix ) {
7784 printf(
"Phase ID= %d \n", npix);
7785 printf(
"Number accessible from first surface = %ld \n", ntop);
7786 printf(
"Number contained in through pathways= %ld \n", nthrough);
7794 alpha_burn = 1. - ( mass_burn /
cemmass );
7801 if ( nthrough > 0 ) {
7813 ( void ) alpha_burn;
7824int CemhydMatStatus :: burnset(
int d1,
int d2,
int d3)
7826 long int ntop, nthrough, icur, inew, ncur, nnew, ntot, count_solid;
7827 int i, j, k, setyet;
7829 int *nmatx, *nmaty, *nmatz;
7830 int xl, xh, j1, k1, px, py, pz, qx, qy, qz;
7831 int xcn, ycn, zcn, x1, y1, z1, igood;
7833 int *nnewx, *nnewy, *nnewz;
7835 float mass_burn = 0.0, alpha_burn = 0.0, con_frac;
7852 for ( k = 0; k <
SYSIZE; k++ ) {
7853 for ( j = 0; j <
SYSIZE; j++ ) {
7854 for ( i = 0; i <
SYSIZE; i++ ) {
7855 newmat [ i ] [ j ] [ k ] =
mic [ i ] [ j ] [ k ];
7865 for ( k = 0; k <
SYSIZE; k++ ) {
7866 for ( j = 0; j <
SYSIZE; j++ ) {
7871 px =
cx(i, j, k, d1, d2, d3);
7872 py =
cy(i, j, k, d1, d2, d3);
7873 pz =
cz(i, j, k, d1, d2, d3);
7876 if ( (
mic [ px ] [ py ] [ pz ] ==
C3S ) ||
7877 (
mic [ px ] [ py ] [ pz ] ==
C2S ) ||
7878 (
mic [ px ] [ py ] [ pz ] ==
SLAG ) ||
7879 (
mic [ px ] [ py ] [ pz ] ==
ASG ) ||
7880 (
mic [ px ] [ py ] [ pz ] ==
CAS2 ) ||
7881 (
mic [ px ] [ py ] [ pz ] ==
POZZ ) ||
7882 (
mic [ px ] [ py ] [ pz ] ==
CSH ) ||
7883 (
mic [ px ] [ py ] [ pz ] ==
C3AH6 ) ||
7884 (
mic [ px ] [ py ] [ pz ] ==
ETTR ) ||
7886 (
mic [ px ] [ py ] [ pz ] ==
C3A ) ||
7887 (
mic [ px ] [ py ] [ pz ] ==
C4AF ) ) {
7900 for ( inew = 1; inew <= ncur; inew++ ) {
7901 xcn = nmatx [ inew ];
7902 ycn = nmaty [ inew ];
7903 zcn = nmatz [ inew ];
7905 qx =
cx(xcn, ycn, zcn, d1, d2, d3);
7906 qy =
cy(xcn, ycn, zcn, d1, d2, d3);
7907 qz =
cz(xcn, ycn, zcn, d1, d2, d3);
7910 for ( jnew = 1; jnew <= 6; jnew++ ) {
7941 }
else if ( y1 < 0 ) {
7948 }
else if ( z1 < 0 ) {
7953 if ( ( x1 >= 0 ) && ( x1 <
SYSIZE ) ) {
7954 px =
cx(x1, y1, z1, d1, d2, d3);
7955 py =
cy(x1, y1, z1, d1, d2, d3);
7956 pz =
cz(x1, y1, z1, d1, d2, d3);
7959 if ( (
mic [ px ] [ py ] [ pz ] ==
CSH ) || (
mic [ px ] [ py ] [ pz ] ==
ETTRC4AF ) || (
mic [ px ] [ py ] [ pz ] ==
C3AH6 ) || (
mic [ px ] [ py ] [ pz ] ==
ETTR ) ) {
7964 printf(
"error in size of nnew %ld\n", nnew);
7967 nnewx [ nnew ] = x1;
7968 nnewy [ nnew ] = y1;
7969 nnewz [ nnew ] = z1;
7972 else if ( ( ( newmat [ qx ] [ qy ] [ qz ] ==
CSH ) || ( newmat [ qx ] [ qy ] [ qz ] ==
ETTRC4AF ) || ( newmat [ qx ] [ qy ] [ qz ] ==
C3AH6 ) || ( newmat [ qx ] [ qy ] [ qz ] ==
ETTR ) ) &&
7973 ( (
mic [ px ] [ py ] [ pz ] ==
C3S ) ||
7974 (
mic [ px ] [ py ] [ pz ] ==
C2S ) ||
7975 (
mic [ px ] [ py ] [ pz ] ==
CAS2 ) ||
7976 (
mic [ px ] [ py ] [ pz ] ==
SLAG ) ||
7977 (
mic [ px ] [ py ] [ pz ] ==
POZZ ) ||
7978 (
mic [ px ] [ py ] [ pz ] ==
ASG ) ||
7979 (
mic [ px ] [ py ] [ pz ] ==
C3A ) ||
7980 (
mic [ px ] [ py ] [ pz ] ==
C4AF ) ) ) {
7985 printf(
"error in size of nnew %ld\n", nnew);
7988 nnewx [ nnew ] = x1;
7989 nnewy [ nnew ] = y1;
7990 nnewz [ nnew ] = z1;
7995 else if ( (
micpart [ qx ] [ qy ] [ qz ] ==
micpart [ px ] [ py ] [ pz ] ) &&
7996 ( (
mic [ px ] [ py ] [ pz ] ==
C3S ) ||
7997 (
mic [ px ] [ py ] [ pz ] ==
C2S ) ||
7998 (
mic [ px ] [ py ] [ pz ] ==
POZZ ) ||
7999 (
mic [ px ] [ py ] [ pz ] ==
SLAG ) ||
8000 (
mic [ px ] [ py ] [ pz ] ==
ASG ) ||
8001 (
mic [ px ] [ py ] [ pz ] ==
CAS2 ) ||
8002 (
mic [ px ] [ py ] [ pz ] ==
C3A ) ||
8003 (
mic [ px ] [ py ] [ pz ] ==
C4AF ) ) && ( ( newmat [ qx ] [ qy ] [ qz ] ==
C3S ) ||
8004 ( newmat [ qx ] [ qy ] [ qz ] ==
C2S ) ||
8005 ( newmat [ qx ] [ qy ] [ qz ] ==
SLAG ) ||
8006 ( newmat [ qx ] [ qy ] [ qz ] ==
ASG ) ||
8007 ( newmat [ qx ] [ qy ] [ qz ] ==
POZZ ) ||
8008 ( newmat [ qx ] [ qy ] [ qz ] ==
CAS2 ) ||
8009 ( newmat [ qx ] [ qy ] [ qz ] ==
C3A ) ||
8010 ( newmat [ qx ] [ qy ] [ qz ] ==
C4AF ) ) ) {
8015 printf(
"error in size of nnew %ld\n", nnew);
8018 nnewx [ nnew ] = x1;
8019 nnewy [ nnew ] = y1;
8020 nnewz [ nnew ] = z1;
8034 for ( icur = 1; icur <= ncur; icur++ ) {
8035 nmatx [ icur ] = nnewx [ icur ];
8036 nmaty [ icur ] = nnewy [ icur ];
8037 nmatz [ icur ] = nnewz [ icur ];
8040 }
while ( nnew > 0 );
8046 for ( j1 = 0; j1 <
SYSIZE; j1++ ) {
8047 for ( k1 = 0; k1 <
SYSIZE; k1++ ) {
8048 px =
cx(xl, j1, k1, d1, d2, d3);
8049 py =
cy(xl, j1, k1, d1, d2, d3);
8050 pz =
cz(xl, j1, k1, d1, d2, d3);
8051 qx =
cx(xh, j1, k1, d1, d2, d3);
8052 qy =
cy(xh, j1, k1, d1, d2, d3);
8053 qz =
cz(xh, j1, k1, d1, d2, d3);
8054 if ( (
mic [ px ] [ py ] [ pz ] ==
BURNT ) && (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) ) {
8058 if (
mic [ px ] [ py ] [ pz ] ==
BURNT ) {
8059 mic [ px ] [ py ] [ pz ] =
BURNT + 1;
8062 if (
mic [ qx ] [ qy ] [ qz ] ==
BURNT ) {
8063 mic [ qx ] [ qy ] [ qz ] =
BURNT + 1;
8076 printf(
"Phase ID= Solid Phases \n");
8077 printf(
"Number accessible from first surface = %ld \n", ntop);
8078 printf(
"Number contained in through pathways= %ld \n", nthrough);
8085 alpha_burn = 1. - ( mass_burn /
cemmass );
8090 if ( count_solid > 0 ) {
8091 con_frac = ( float ) nthrough / (
float ) count_solid;
8099 if ( con_frac > 0.985 ) {
8106 for ( i = 0; i <
SYSIZE; i++ ) {
8107 for ( j = 0; j <
SYSIZE; j++ ) {
8108 for ( k = 0; k <
SYSIZE; k++ ) {
8109 if (
mic [ i ] [ j ] [ k ] >=
BURNT ) {
8110 mic [ i ] [ j ] [ k ] = newmat [ i ] [ j ] [ k ];
8128 ( void ) alpha_burn;
8134void CemhydMatStatus :: parthyd()
8136 int norig [ 100000 ], nleft [ 100000 ];
8138 char valmic, valmicorig;
8139 int valpart, partmax;
8144 for ( ix = 0; ix < 100000; ix++ ) {
8145 nleft [ ix ] = norig [ ix ] = 0;
8148 phydfile = fopen(
phrname,
"a");
8153 for ( ix = 0; ix <
SYSIZE; ix++ ) {
8155 for ( iz = 0; iz <
SYSIZE; iz++ ) {
8156 if (
micpart [ ix ] [
iy ] [ iz ] != 0 ) {
8158 if ( valpart > partmax ) {
8162 valmic =
mic [ ix ] [
iy ] [ iz ];
8163 if ( ( valmic ==
C3S ) || ( valmic ==
C2S ) || ( valmic ==
C3A ) || ( valmic ==
C4AF ) ) {
8164 nleft [ valpart ] += 1;
8167 valmicorig =
micorig [ ix ] [
iy ] [ iz ];
8168 if ( ( valmicorig ==
C3S ) || ( valmicorig ==
C2S ) || ( valmicorig ==
C3A ) || ( valmicorig ==
C4AF ) ) {
8169 norig [ valpart ] += 1;
8177 for ( ix = 100; ix <= partmax; ix++ ) {
8179 if ( norig [ ix ] != 0 ) {
8180 alpart = 1. - ( float ) nleft [ ix ] / (
float ) norig [ ix ];
8183 fprintf(phydfile,
"%d %d %d %.3f\n", ix, norig [ ix ], nleft [ ix ], alpart);
8192int CemhydMatStatus :: moveone(
int *xloc,
int *yloc,
int *zloc,
int *act,
int sumold)
8194 int plok, sumnew, xl1, yl1, zl1, act1;
8205 plok = ( int ) ( 6. *
ran1(
seed) );
8206 if ( ( plok > 5 ) || ( plok < 0 ) ) {
8218 if ( sumold % 2 != 0 ) {
8230 if ( sumold % 3 != 0 ) {
8242 if ( sumold % 5 != 0 ) {
8254 if ( sumold % 7 != 0 ) {
8266 if ( sumold % 11 != 0 ) {
8278 if ( sumold % 13 != 0 ) {
8303int CemhydMatStatus :: edgecnt(
int xck,
int yck,
int zck,
int ph1,
int ph2,
int ph3)
8305 int ixe, iye, ize, edgeback, x2, y2, z2, check;
8312 for ( ixe = ( -1 ); ixe <= 1; ixe++ ) {
8314 for ( iye = ( -1 ); iye <= 1; iye++ ) {
8316 for ( ize = ( -1 ); ize <= 1; ize++ ) {
8317 if ( ( ixe != 0 ) || ( iye != 0 ) || ( ize != 0 ) ) {
8322 }
else if ( x2 >=
SYSIZE ) {
8328 }
else if ( y2 >=
SYSIZE ) {
8334 }
else if ( z2 >=
SYSIZE ) {
8338 check =
mic [ x2 ] [ y2 ] [ z2 ];
8339 if ( ( check != ph1 ) && ( check != ph2 ) && ( check != ph3 ) ) {
8348 return ( edgeback );
8354void CemhydMatStatus :: extcsh()
8356 int numnear, xchr, ychr, zchr, fchr, check, msface;
8363 while ( fchr == 0 ) {
8381 check =
mic [ xchr ] [ ychr ] [ zchr ];
8388 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8389 mic [ xchr ] [ ychr ] [ zchr ] =
CSH;
8394 msface = ( int ) ( 3. *
ran1(
seed) + 1. );
8399 faces [ xchr ] [ ychr ] [ zchr ] = msface;
8415int CemhydMatStatus :: movecsh(
int xcur,
int ycur,
int zcur,
int finalstep,
int cycorig)
8417 int xnew, ynew, znew, action, sumback, sumin, check;
8418 int msface, mstest = 0, mstest2 = 0;
8419 float prcsh, prcsh1, prcsh2, prtest;
8427 sumback =
moveone(& xnew, & ynew, & znew, & action, sumin);
8430 if ( xnew != xcur ) {
8435 if ( ynew != ycur ) {
8440 if ( znew != zcur ) {
8446 if ( action == 0 ) {
8447 printf(
"Error in value of action \n");
8450 check =
mic [ xnew ] [ ynew ] [ znew ];
8456 if ( ( check ==
CSH ) && ( (
cshgeom == 0 ) || (
faces [ xnew ] [ ynew ] [ znew ] == 0 ) || (
faces [ xnew ] [ ynew ] [ znew ] == mstest ) || (
faces [ xnew ] [ ynew ] [ znew ] == mstest2 ) ) ) {
8462 if ( prcsh1 <= prtest ) {
8463 mic [ xcur ] [ ycur ] [ zcur ] =
CSH;
8465 faces [ xcur ] [ ycur ] [ zcur ] =
faces [ xnew ] [ ynew ] [ znew ];
8477 if ( prtest > 1.0 ) {
8479 if ( prcsh2 < ( prtest - 1.0 ) ) {
8487 else if ( ( check ==
SLAGCSH ) || ( check ==
POZZCSH ) || ( finalstep == 1 ) ||
8488 ( ( ( check ==
C3S ) || ( check ==
C2S ) ) && ( prcsh < 0.001 ) ) ||
8489 ( ( ( check ==
C3A ) || ( check ==
C4AF ) ) && ( prcsh < 0.2 ) ) ||
8490 ( ( check ==
CH ) && ( prcsh < 0.01 ) ) ||
8497 if ( prcsh1 <= prtest ) {
8498 mic [ xcur ] [ ycur ] [ zcur ] =
CSH;
8501 msface = ( int ) ( 2. *
ran1(
seed) + 1. );
8506 if ( msface == 1 ) {
8507 faces [ xcur ] [ ycur ] [ zcur ] = mstest;
8509 faces [ xcur ] [ ycur ] [ zcur ] = mstest2;
8522 if ( prtest > 1.0 ) {
8524 if ( prcsh2 < ( prtest - 1.0 ) ) {
8532 if ( action != 0 ) {
8553void CemhydMatStatus :: extfh3(
int xpres,
int ypres,
int zpres)
8555 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
8564 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
8570 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
8571 if ( newact == 0 ) {
8572 printf(
"Error in value of newact in extfh3 \n");
8575 check =
mic [ xchr ] [ ychr ] [ zchr ];
8580 mic [ xchr ] [ ychr ] [ zchr ] =
FH3;
8592 while ( fchr == 0 ) {
8610 check =
mic [ xchr ] [ ychr ] [ zchr ];
8616 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8617 mic [ xchr ] [ ychr ] [ zchr ] =
FH3;
8634int CemhydMatStatus :: extettr(
int xpres,
int ypres,
int zpres,
int etype)
8636 int check, newact, multf, numnear, sump, xchr, ychr, zchr, fchr, i1;
8637 int numalum, numsil;
8638 float pneigh, ptest;
8648 for ( i1 = 1; ( ( i1 <= 1000 ) && ( fchr == 0 ) ); i1++ ) {
8654 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
8655 if ( newact == 0 ) {
8656 printf(
"Error in value of action \n");
8659 check =
mic [ xchr ] [ ychr ] [ zchr ];
8667 numsil = 26 - numsil;
8671 numalum = 26 - numalum;
8675 numalum = 26 - numalum;
8678 pneigh = ( float ) ( numnear + 1 ) / 26.0;
8680 if ( numalum <= 1 ) {
8684 if ( numalum >= 2 ) {
8688 if ( numalum >= 3 ) {
8692 if ( numalum >= 5 ) {
8698 if ( pneigh >= ptest ) {
8700 mic [ xchr ] [ ychr ] [ zchr ] =
ETTR;
8718 while ( fchr == 0 ) {
8737 check =
mic [ xchr ] [ ychr ] [ zchr ];
8741 numsil = 26 - numsil;
8750 if ( ( tries > 5000 ) || ( ( numnear < 26 ) && ( numsil < 1 ) ) ) {
8752 mic [ xchr ] [ ychr ] [ zchr ] =
ETTR;
8774void CemhydMatStatus :: extch()
8776 int numnear, xchr, ychr, zchr, fchr, check;
8783 while ( fchr == 0 ) {
8801 check =
mic [ xchr ] [ ychr ] [ zchr ];
8808 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8809 mic [ xchr ] [ ychr ] [ zchr ] =
CH;
8821void CemhydMatStatus :: extgyps(
int xpres,
int ypres,
int zpres)
8823 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
8832 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
8838 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
8839 if ( newact == 0 ) {
8840 printf(
"Error in value of newact in extfh3 \n");
8843 check =
mic [ xchr ] [ ychr ] [ zchr ];
8860 while ( fchr == 0 ) {
8878 check =
mic [ xchr ] [ ychr ] [ zchr ];
8884 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
8900int CemhydMatStatus :: moveanh(
int xcur,
int ycur,
int zcur,
int finalstep,
float nucprgyp)
8902 int xnew, ynew, znew, action, sumback, sumin, check = -1;
8903 int nexp, iexp, xexp, yexp, zexp, newact, ettrtype;
8904 float pgen, pexp, pext, p2diff;
8910 if ( ( nucprgyp >= pgen ) || ( finalstep == 1 ) ) {
8925 sumback =
moveone(& xnew, & ynew, & znew, & action, sumin);
8927 if ( action == 0 ) {
8928 printf(
"Error in value of action \n");
8931 check =
mic [ xnew ] [ ynew ] [ znew ];
8953 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
8961 count [ check ] -= 1;
8968 if ( pexp <= 0.569 ) {
8969 if ( ettrtype == 0 ) {
8970 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
8981 if ( check ==
C3A ) {
8982 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
8985 if ( ettrtype == 0 ) {
9003 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9004 newact =
extettr(xexp, yexp, zexp, ettrtype);
9056 if ( pexp <= 0.6935 ) {
9057 newact =
extettr(xexp, yexp, zexp, ettrtype);
9073 if ( pexp <= 0.8174 ) {
9080 if ( pext < 0.2584 ) {
9086 if ( pext < 0.5453 ) {
9087 extfh3(xnew, ynew, znew);
9092 mic [ xnew ] [ ynew ] [ znew ] =
C4AF;
9102 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9103 newact =
extettr(xexp, yexp, zexp, 1);
9155 if ( pexp <= 0.6935 ) {
9156 newact =
extettr(xexp, yexp, zexp, 1);
9163 if ( action != 0 ) {
9186int CemhydMatStatus :: movehem(
int xcur,
int ycur,
int zcur,
int finalstep,
float nucprgyp)
9188 int xnew, ynew, znew, action, sumback, sumin, check = -1;
9189 int nexp, iexp, xexp, yexp, zexp, newact, ettrtype;
9190 float pgen, pexp, pext, p2diff;
9196 if ( ( nucprgyp >= pgen ) || ( finalstep == 1 ) ) {
9212 sumback =
moveone(& xnew, & ynew, & znew, & action, sumin);
9214 if ( action == 0 ) {
9215 printf(
"Error in value of action \n");
9218 check =
mic [ xnew ] [ ynew ] [ znew ];
9240 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
9248 count [ check ] -= 1;
9255 if ( pexp <= 0.5583 ) {
9256 if ( ettrtype == 0 ) {
9257 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
9268 if ( check ==
C3A ) {
9269 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
9272 if ( ettrtype == 0 ) {
9290 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9291 newact =
extettr(xexp, yexp, zexp, ettrtype);
9343 if ( pexp <= 0.6053 ) {
9344 newact =
extettr(xexp, yexp, zexp, ettrtype);
9360 if ( pexp <= 0.802 ) {
9367 if ( pext < 0.2584 ) {
9373 if ( pext < 0.5453 ) {
9374 extfh3(xnew, ynew, znew);
9379 mic [ xnew ] [ ynew ] [ znew ] =
C4AF;
9389 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9390 newact =
extettr(xexp, yexp, zexp, 1);
9442 if ( pexp <= 0.6053 ) {
9443 newact =
extettr(xexp, yexp, zexp, 1);
9450 if ( action != 0 ) {
9471int CemhydMatStatus :: extfreidel(
int xpres,
int ypres,
int zpres)
9473 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
9482 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
9488 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
9489 if ( newact == 0 ) {
9490 printf(
"Error in value of newact in extfreidel \n");
9493 check =
mic [ xchr ] [ ychr ] [ zchr ];
9510 while ( fchr == 0 ) {
9529 check =
mic [ xchr ] [ ychr ] [ zchr ];
9535 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
9552int CemhydMatStatus :: extstrat(
int xpres,
int ypres,
int zpres)
9554 int multf, numnear, sump, xchr, ychr, zchr, check, fchr, i1, newact;
9563 for ( i1 = 1; ( ( i1 <= 500 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
9569 multf =
moveone(& xchr, & ychr, & zchr, & newact, sump);
9570 if ( newact == 0 ) {
9571 printf(
"Error in value of newact in extstrat \n");
9574 check =
mic [ xchr ] [ ychr ] [ zchr ];
9579 mic [ xchr ] [ ychr ] [ zchr ] =
STRAT;
9591 while ( fchr == 0 ) {
9610 check =
mic [ xchr ] [ ychr ] [ zchr ];
9616 if ( ( numnear < 26 ) || ( tries > 5000 ) ) {
9617 mic [ xchr ] [ ychr ] [ zchr ] =
STRAT;
9633int CemhydMatStatus :: movegyp(
int xcur,
int ycur,
int zcur,
int finalstep)
9635 int check, xnew, ynew, znew, action, nexp, iexp;
9636 int xexp, yexp, zexp, newact, sumold, sumgarb, ettrtype;
9637 float pexp, pext, p2diff;
9643 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFGYP ) {
9653 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
9654 if ( action == 0 ) {
9655 printf(
"Error in value of action in movegyp \n");
9658 check =
mic [ xnew ] [ ynew ] [ znew ];
9667 mic [ xcur ] [ ycur ] [ zcur ] =
ABSGYP;
9678 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
9686 count [ check ] -= 1;
9693 if ( pexp <= 0.40 ) {
9694 if ( ettrtype == 0 ) {
9695 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
9706 if ( check ==
C3A ) {
9707 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
9710 if ( ettrtype == 0 ) {
9728 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9729 newact =
extettr(xexp, yexp, zexp, ettrtype);
9781 if ( pexp <= 0.30 ) {
9782 newact =
extettr(xexp, yexp, zexp, ettrtype);
9798 if ( pexp <= 0.575 ) {
9805 if ( pext < 0.2584 ) {
9811 if ( pext < 0.5453 ) {
9812 extfh3(xnew, ynew, znew);
9817 mic [ xnew ] [ ynew ] [ znew ] =
C4AF;
9827 for ( iexp = 1; iexp <= nexp; iexp++ ) {
9828 newact =
extettr(xexp, yexp, zexp, 1);
9880 if ( pexp <= 0.30 ) {
9881 newact =
extettr(xexp, yexp, zexp, 1);
9889 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
9893 mic [ xcur ] [ ycur ] [ zcur ] =
GYPSUM;
9896 if ( action != 0 ) {
9918int CemhydMatStatus :: movecacl2(
int xcur,
int ycur,
int zcur,
int finalstep)
9920 int check, xnew, ynew, znew, action, nexp, iexp;
9921 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
9939 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
9940 if ( action == 0 ) {
9941 printf(
"Error in value of action in movecacl2 \n");
9944 check =
mic [ xnew ] [ ynew ] [ znew ];
9953 count [ check ] -= 1;
9960 if ( pexp <= 0.5793 ) {
9975 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10028 if ( pexp <= 0.3295 ) {
10034 else if ( check ==
C4AF ) {
10046 if ( pexp <= 0.4033 ) {
10053 if ( pext < 0.6412 ) {
10059 if ( pext < 0.3522 ) {
10060 extfh3(xnew, ynew, znew);
10063 extfh3(xnew, ynew, znew);
10074 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10077 switch ( newact ) {
10127 if ( pexp <= 0.3176 ) {
10136 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10140 mic [ xcur ] [ ycur ] [ zcur ] =
CACL2;
10143 if ( action != 0 ) {
10169int CemhydMatStatus :: movecas2(
int xcur,
int ycur,
int zcur,
int finalstep)
10171 int check, xnew, ynew, znew, action, nexp, iexp;
10172 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
10180 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFCAS2 ) {
10190 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10191 if ( action == 0 ) {
10192 printf(
"Error in value of action in movecas2 \n");
10195 check =
mic [ xnew ] [ ynew ] [ znew ];
10202 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
10211 if ( pexp <= 0.886 ) {
10212 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
10214 count [ check ] -= 1;
10223 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10224 newact =
extstrat(xexp, yexp, zexp);
10226 switch ( newact ) {
10276 if ( pexp <= 0.286 ) {
10277 newact =
extstrat(xexp, yexp, zexp);
10282 else if ( check ==
C4AF ) {
10283 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
10294 if ( pexp <= 0.786 ) {
10295 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
10302 if ( pext < 0.329 ) {
10309 if ( pext < 0.6938 ) {
10310 extfh3(xnew, ynew, znew);
10322 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10323 newact =
extstrat(xexp, yexp, zexp);
10325 switch ( newact ) {
10375 if ( pexp <= 0.37 ) {
10376 newact =
extstrat(xexp, yexp, zexp);
10384 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10388 mic [ xcur ] [ ycur ] [ zcur ] =
CAS2;
10391 if ( action != 0 ) {
10417int CemhydMatStatus :: moveas(
int xcur,
int ycur,
int zcur,
int finalstep)
10419 int check, xnew, ynew, znew, action, nexp, iexp;
10420 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
10428 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFAS ) {
10438 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10439 if ( action == 0 ) {
10440 printf(
"Error in value of action in moveas \n");
10443 check =
mic [ xnew ] [ ynew ] [ znew ];
10447 if ( ( check ==
CH ) || ( check ==
DIFFCH ) ) {
10450 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
10452 count [ check ] -= 1;
10459 if ( pexp <= 0.7538 ) {
10460 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
10474 for ( iexp = 1; iexp <= nexp; iexp++ ) {
10475 newact =
extstrat(xexp, yexp, zexp);
10477 switch ( newact ) {
10527 if ( pexp <= 0.326 ) {
10528 newact =
extstrat(xexp, yexp, zexp);
10534 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10538 mic [ xcur ] [ ycur ] [ zcur ] =
ASG;
10541 if ( action != 0 ) {
10545 mic [ xnew ] [ ynew ] [ znew ] =
DIFFAS;
10567int CemhydMatStatus :: movecaco3(
int xcur,
int ycur,
int zcur,
int finalstep)
10569 int check, xnew, ynew, znew, action;
10570 int xexp, yexp, zexp, newact, sumold, sumgarb, keep;
10588 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10589 if ( action == 0 ) {
10590 printf(
"Error in value of action in moveas \n");
10593 check =
mic [ xnew ] [ ynew ] [ znew ];
10598 if ( check ==
AFM ) {
10602 if ( pexp <= 0.479192 ) {
10603 mic [ xnew ] [ ynew ] [ znew ] =
AFMC;
10606 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
10610 count [ check ] -= 1;
10616 if ( pexp <= 0.078658 ) {
10617 mic [ xcur ] [ ycur ] [ zcur ] =
AFMC;
10632 if ( pexp <= 0.26194 ) {
10633 newact =
extettr(xexp, yexp, zexp, 0);
10639 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10643 mic [ xcur ] [ ycur ] [ zcur ] =
CACO3;
10646 if ( action != 0 ) {
10672void CemhydMatStatus :: extafm(
int xpres,
int ypres,
int zpres)
10674 int check, sump, xchr, ychr, zchr, fchr, i1, newact, numnear;
10683 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
10689 sump *=
moveone(& xchr, & ychr, & zchr, & newact, sump);
10690 if ( newact == 0 ) {
10691 printf(
"Error in value of newact in extafm \n");
10694 check =
mic [ xchr ] [ ychr ] [ zchr ];
10698 mic [ xchr ] [ ychr ] [ zchr ] =
AFM;
10708 while ( fchr == 0 ) {
10726 check =
mic [ xchr ] [ ychr ] [ zchr ];
10733 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
10734 mic [ xchr ] [ ychr ] [ zchr ] =
AFM;
10747int CemhydMatStatus :: moveettr(
int xcur,
int ycur,
int zcur,
int finalstep)
10749 int check, xnew, ynew, znew, action;
10750 int sumold, sumgarb;
10751 float pexp, pafm, pgrow;
10755 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFETTR ) {
10766 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10767 if ( action == 0 ) {
10768 printf(
"Error in value of action in moveettr \n");
10771 check =
mic [ xnew ] [ ynew ] [ znew ];
10775 if ( check ==
C4AF ) {
10777 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
10787 if ( pexp <= 0.278 ) {
10788 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
10793 if ( pafm <= 0.3241 ) {
10799 if ( pafm <= 0.4313 ) {
10800 extfh3(xnew, ynew, znew);
10802 }
else if ( pexp <= 0.348 ) {
10803 mic [ xnew ] [ ynew ] [ znew ] =
FH3;
10812 else if ( ( check ==
C3A ) || ( check ==
DIFFC3A ) ) {
10815 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
10818 count [ check ] -= 1;
10824 if ( pexp <= 0.2424 ) {
10825 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
10831 if ( check ==
C3A ) {
10832 mic [ xnew ] [ ynew ] [ znew ] =
C3A;
10845 if ( pexp <= pafm ) {
10846 extafm(xcur, ycur, zcur);
10850 else if ( check ==
ETTR ) {
10853 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
10862 if ( ( action != 0 ) && ( finalstep == 1 ) ) {
10866 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
10869 if ( action != 0 ) {
10890void CemhydMatStatus :: extpozz(
int xpres,
int ypres,
int zpres)
10892 int check, sump, xchr, ychr, zchr, fchr, i1, action, numnear;
10901 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
10907 sump *=
moveone(& xchr, & ychr, & zchr, & action, sump);
10908 if ( action == 0 ) {
10909 printf(
"Error in value of action in extpozz \n");
10912 check =
mic [ xchr ] [ ychr ] [ zchr ];
10926 while ( fchr == 0 ) {
10944 check =
mic [ xchr ] [ ychr ] [ zchr ];
10950 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
10964int CemhydMatStatus :: movefh3(
int xcur,
int ycur,
int zcur,
int finalstep,
float nucprob)
10966 int check, xnew, ynew, znew, action, sumold, sumgarb;
10972 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
10974 mic [ xcur ] [ ycur ] [ zcur ] =
FH3;
10984 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
10985 if ( action == 0 ) {
10986 printf(
"Error in value of action in movefh3 \n");
10989 check =
mic [ xnew ] [ ynew ] [ znew ];
10992 if ( check ==
FH3 ) {
10993 mic [ xcur ] [ ycur ] [ zcur ] =
FH3;
10999 if ( action != 0 ) {
11021int CemhydMatStatus :: movech(
int xcur,
int ycur,
int zcur,
int finalstep,
float nucprob)
11023 int check, xnew, ynew, znew, action, sumgarb, sumold;
11024 float pexp, pgen, pfix;
11028 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
11030 mic [ xcur ] [ ycur ] [ zcur ] =
CH;
11040 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
11041 if ( action == 0 ) {
11042 printf(
"Error in value of action in movech \n");
11045 check =
mic [ xnew ] [ ynew ] [ znew ];
11048 if ( ( check ==
CH ) && ( pgen <=
CHGROW ) ) {
11049 mic [ xcur ] [ ycur ] [ zcur ] =
CH;
11057 mic [ xcur ] [ ycur ] [ zcur ] =
CH;
11064 else if ( ( pgen <=
ppozz ) && ( check ==
POZZ ) && (
npr <= (
int ) ( (
float )
nfill * 1.35 ) ) ) {
11074 if ( pfix <= ( 1. / 1.35 ) ) {
11086 if ( pexp <= 0.05466 ) {
11089 }
else if ( check ==
DIFFAS ) {
11091 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
11099 if ( pfix <= 0.7538 ) {
11100 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
11109 if ( pexp <= 0.5035 ) {
11114 if ( action != 0 ) {
11118 mic [ xnew ] [ ynew ] [ znew ] =
DIFFCH;
11136void CemhydMatStatus :: extc3ah6(
int xpres,
int ypres,
int zpres)
11138 int check, sump, xchr, ychr, zchr, fchr, i1, action, numnear;
11147 for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
11153 sump *=
moveone(& xchr, & ychr, & zchr, & action, sump);
11154 if ( action == 0 ) {
11155 printf(
"Error in action value in extc3ah6 \n");
11158 check =
mic [ xchr ] [ ychr ] [ zchr ];
11162 mic [ xchr ] [ ychr ] [ zchr ] =
C3AH6;
11171 while ( fchr == 0 ) {
11188 check =
mic [ xchr ] [ ychr ] [ zchr ];
11194 if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
11195 mic [ xchr ] [ ychr ] [ zchr ] =
C3AH6;
11208int CemhydMatStatus :: movec3a(
int xcur,
int ycur,
int zcur,
int finalstep,
float nucprob)
11210 int check, xnew, ynew, znew, action, sumgarb, sumold;
11211 int xexp, yexp, zexp, nexp, iexp, newact;
11212 float pgen, pexp, pafm, pgrow, p2diff;
11215 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFC3A ) {
11224 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
11226 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11233 if ( pexp <= 0.69 ) {
11243 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
11244 if ( action == 0 ) {
11245 printf(
"Error in value of action in movec3a \n");
11248 check =
mic [ xnew ] [ ynew ] [ znew ];
11251 if ( check ==
C3AH6 ) {
11256 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11263 if ( pexp <= 0.69 ) {
11272 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
11281 if ( pexp <= 0.40 ) {
11282 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
11298 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11299 newact =
extettr(xexp, yexp, zexp, 0);
11301 switch ( newact ) {
11351 if ( pexp <= 0.30 ) {
11352 newact =
extettr(xexp, yexp, zexp, 0);
11359 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
11368 if ( pexp <= 0.5583 ) {
11369 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
11385 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11386 newact =
extettr(xexp, yexp, zexp, 0);
11388 switch ( newact ) {
11438 if ( pexp <= 0.6053 ) {
11439 newact =
extettr(xexp, yexp, zexp, 0);
11446 mic [ xnew ] [ ynew ] [ znew ] =
ETTR;
11455 if ( pexp <= 0.569 ) {
11456 mic [ xcur ] [ ycur ] [ zcur ] =
ETTR;
11472 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11473 newact =
extettr(xexp, yexp, zexp, 0);
11475 switch ( newact ) {
11525 if ( pexp <= 0.6935 ) {
11526 newact =
extettr(xexp, yexp, zexp, 0);
11542 if ( pexp <= 0.5793 ) {
11557 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11560 switch ( newact ) {
11610 if ( pexp <= 0.3295 ) {
11618 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
11627 if ( pexp <= 0.886 ) {
11628 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
11643 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11644 newact =
extstrat(xexp, yexp, zexp);
11646 switch ( newact ) {
11696 if ( pexp <= 0.286 ) {
11697 newact =
extstrat(xexp, yexp, zexp);
11708 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
11711 count [ check ] -= 1;
11716 if ( pexp <= 0.2424 ) {
11717 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
11728 if ( pexp <= pafm ) {
11729 extafm(xnew, ynew, znew);
11733 if ( ( action != 0 ) && ( action != 7 ) ) {
11755int CemhydMatStatus :: movec4a(
int xcur,
int ycur,
int zcur,
int finalstep,
float nucprob)
11757 int check, xnew, ynew, znew, action, sumgarb, sumold;
11758 int xexp, yexp, zexp, nexp, iexp, newact;
11759 float pgen, pexp, pafm, pgrow, p2diff;
11762 if (
mic [ xcur ] [ ycur ] [ zcur ] !=
DIFFC4A ) {
11771 if ( ( nucprob >= pgen ) || ( finalstep == 1 ) ) {
11773 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11780 if ( pexp <= 0.69 ) {
11790 sumgarb =
moveone(& xnew, & ynew, & znew, & action, sumold);
11791 if ( action == 0 ) {
11792 printf(
"Error in value of action in movec4a \n");
11795 check =
mic [ xnew ] [ ynew ] [ znew ];
11798 if ( check ==
C3AH6 ) {
11803 mic [ xcur ] [ ycur ] [ zcur ] =
C3AH6;
11810 if ( pexp <= 0.69 ) {
11828 if ( pexp <= 0.40 ) {
11845 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11846 newact =
extettr(xexp, yexp, zexp, 1);
11848 switch ( newact ) {
11898 if ( pexp <= 0.30 ) {
11899 newact =
extettr(xexp, yexp, zexp, 1);
11915 if ( pexp <= 0.5583 ) {
11932 for ( iexp = 1; iexp <= nexp; iexp++ ) {
11933 newact =
extettr(xexp, yexp, zexp, 1);
11935 switch ( newact ) {
11985 if ( pexp <= 0.6053 ) {
11986 newact =
extettr(xexp, yexp, zexp, 1);
12002 if ( pexp <= 0.569 ) {
12019 for ( iexp = 1; iexp <= nexp; iexp++ ) {
12020 newact =
extettr(xexp, yexp, zexp, 1);
12022 switch ( newact ) {
12072 if ( pexp <= 0.6935 ) {
12073 newact =
extettr(xexp, yexp, zexp, 1);
12089 if ( pexp <= 0.5793 ) {
12104 for ( iexp = 1; iexp <= nexp; iexp++ ) {
12107 switch ( newact ) {
12157 if ( pexp <= 0.3295 ) {
12165 mic [ xnew ] [ ynew ] [ znew ] =
STRAT;
12174 if ( pexp <= 0.886 ) {
12175 mic [ xcur ] [ ycur ] [ zcur ] =
STRAT;
12190 for ( iexp = 1; iexp <= nexp; iexp++ ) {
12191 newact =
extstrat(xexp, yexp, zexp);
12193 switch ( newact ) {
12243 if ( pexp <= 0.286 ) {
12244 newact =
extstrat(xexp, yexp, zexp);
12255 mic [ xnew ] [ ynew ] [ znew ] =
AFM;
12258 count [ check ] -= 1;
12263 if ( pexp <= 0.2424 ) {
12264 mic [ xcur ] [ ycur ] [ zcur ] =
AFM;
12275 if ( pexp <= pafm ) {
12276 extafm(xnew, ynew, znew);
12280 if ( ( action != 0 ) && ( action != 7 ) ) {
12302void CemhydMatStatus :: hydrate(
int fincyc,
int stepmax,
float chpar1,
float chpar2,
float hgpar1,
float hgpar2,
float fhpar1,
float fhpar2,
float gypar1,
float gypar2)
12304 int xpl, ypl, zpl, phpl, agepl, xpnew, ypnew, zpnew;
12305 float chprob, c3ah6prob, fh3prob, gypprob;
12306 long int nleft, ntodo, ndale;
12307 int istep, termflag, reactf = -1;
12309 struct ants *curant, *antgone;
12316 for ( istep = 1; ( ( istep <= stepmax ) && ( nleft > 0 ) ); istep++ ) {
12317 if ( ( fincyc == 1 ) && ( istep == stepmax ) ) {
12326 chprob = chpar1 * ( 1. - beterm );
12328 c3ah6prob = hgpar1 * ( 1. - beterm );
12330 fh3prob = fhpar1 * ( 1. - beterm );
12332 gypprob = gypar1 * ( 1. - beterm );
12336 while ( curant != NULL ) {
12342 agepl = curant->cycbirth;
12348 reactf =
movecsh(xpl, ypl, zpl, termflag, agepl);
12349 }
else if ( phpl ==
DIFFANH ) {
12352 reactf =
moveanh(xpl, ypl, zpl, termflag, gypprob);
12353 }
else if ( phpl ==
DIFFHEM ) {
12356 reactf =
movehem(xpl, ypl, zpl, termflag, gypprob);
12357 }
else if ( phpl ==
DIFFCH ) {
12360 reactf =
movech(xpl, ypl, zpl, termflag, chprob);
12361 }
else if ( phpl ==
DIFFFH3 ) {
12364 reactf =
movefh3(xpl, ypl, zpl, termflag, fh3prob);
12365 }
else if ( phpl ==
DIFFGYP ) {
12368 reactf =
movegyp(xpl, ypl, zpl, termflag);
12369 }
else if ( phpl ==
DIFFC3A ) {
12372 reactf =
movec3a(xpl, ypl, zpl, termflag, c3ah6prob);
12373 }
else if ( phpl ==
DIFFC4A ) {
12376 reactf =
movec4a(xpl, ypl, zpl, termflag, c3ah6prob);
12380 reactf =
moveettr(xpl, ypl, zpl, termflag);
12384 reactf =
movecacl2(xpl, ypl, zpl, termflag);
12388 reactf =
movecas2(xpl, ypl, zpl, termflag);
12389 }
else if ( phpl ==
DIFFAS ) {
12392 reactf =
moveas(xpl, ypl, zpl, termflag);
12396 reactf =
movecaco3(xpl, ypl, zpl, termflag);
12398 printf(
"Error in ID of phase \n");
12402 if ( reactf != 0 ) {
12409 switch ( reactf ) {
12419 if ( xpnew >=
SYSIZE ) {
12433 if ( ypnew >=
SYSIZE ) {
12447 if ( zpnew >=
SYSIZE ) {
12461 curant = curant->nextant;
12465 if ( ndale == 1 ) {
12466 headant->nextant = curant->nextant;
12468 ( curant->prevant )->nextant = curant->nextant;
12471 if ( curant->nextant != NULL ) {
12472 ( curant->nextant )->prevant = curant->prevant;
12497 float err, dxold, cdx, abx;
12498 fcomplex_cem sq, h,
gp, gm, g2, g, b, d, dx, f, x1;
12501 for ( iter = 1; iter <=
MAXIT; iter++ ) {
12506 for ( j = m - 1; j >= 0; j-- ) {
12510 err =
Cabs(b) + abx * err;
12514 if (
Cabs(b) <= err ) {
12529 x1 =
Csub(* x, dx);
12530 if ( x->
r == x1.
r && x->
i == x1.
i ) {
12536 if ( iter > 6 && cdx >= dxold ) {
12542 if ( cdx <= eps *
Cabs(* x) ) {
12548 nrerror(
"Too many iterations in routine CemhydMat - LAGUER\n");
12561 for ( j = 0; j <= m; j++ ) {
12562 ad [ j ] = a [ j ];
12565 for ( j = m; j >= 1; j-- ) {
12568 if ( fabs(x.
i) <= ( 2.0 *
EPSP * fabs(x.
r) ) ) {
12574 for ( jj = j - 1; jj >= 0; jj-- ) {
12582 for ( j = 1; j <= m; j++ ) {
12587 for ( j = 2; j <= m; j++ ) {
12589 for ( i = j - 1; i >= 1; i-- ) {
12590 if ( roots [ i ].r <= x.
r ) {
12594 roots [ i + 1 ] = roots [ i ];
12597 roots [ i + 1 ] = x;
12605void CemhydMatStatus :: pHpred()
12607 int j, syngen_change = 0, syn_old = 0, counter = 0;
12608 double concnaplus, conckplus;
12609 double concohminus = 0., A, B, C, conctest, concsulfate1;
12610 double volpore, grams_cement;
12611 double releasedna, releasedk, activitySO4, activityK, test_precip;
12612 double activityCa, activityOH, Istrength, Anow, Bnow, Inew;
12613 double conductivity = 0.0;
12615 float sumbest, sumtest, pozzreact, KspCH;
12635 volpore /= grams_cement;
12642 releasedk /=
MMK2O;
12650 releasedk /=
MMK2O;
12663 activityCa = activityOH = activitySO4 = activityK = 1.0;
12665 if ( ( ( concnaplus + conckplus ) > 0.0 ) && (
soluble [
ETTR ] == 0 ) ) {
12668 if ( Istrength < 1.0 ) {
12672 while ( ( fabs(Istrength - Inew) / Istrength ) > 0.10 ) {
12674 if ( Istrength < 1.0 ) {
12681 activityCa = ( -Anow *
zCa *
zCa * sqrt(Istrength) ) / ( 1. +
aCa * Bnow * sqrt(Istrength) );
12682 activityCa += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zCa *
zCa * Istrength / sqrt(1000.);
12683 activityCa = exp(activityCa);
12684 activityOH = ( -Anow *
zOH *
zOH * sqrt(Istrength) ) / ( 1. +
aOH * Bnow * sqrt(Istrength) );
12685 activityOH += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zOH *
zOH * Istrength / sqrt(1000.);
12686 activityOH = exp(activityOH);
12687 activityK = ( -Anow *
zK *
zK * sqrt(Istrength) ) / ( 1. +
aK * Bnow * sqrt(Istrength) );
12688 activityK += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zK *
zK * Istrength / sqrt(1000.);
12689 activityK = exp(activityK);
12690 activitySO4 = ( -Anow *
zSO4 *
zSO4 * sqrt(Istrength) ) / ( 1. +
aSO4 * Bnow * sqrt(Istrength) );
12691 activitySO4 += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zSO4 *
zSO4 * Istrength / sqrt(1000.);
12692 activitySO4 = exp(activitySO4);
12697 A = ( -KspCH / ( activityCa * activityOH * activityOH ) );
12698 B = conckplus + concnaplus;
12699 C = ( -2. *
KspGypsum / ( activityCa * activitySO4 ) );
12700 concohminus = conckplus + concnaplus;
12715 zroots(coef, 4, roots, 1);
12718 for ( j = 1; j <= 4; j++ ) {
12721 if ( ( ( roots [ j ].i ) == 0.0 ) && ( ( roots [ j ].r ) > 0.0 ) ) {
12722 conctest = sqrt( KspCH / ( roots [ j ].r * activityCa * activityOH * activityOH ) );
12723 concsulfate1 =
KspGypsum / ( roots [ j ].
r * activityCa * activitySO4 );
12724 sumtest = concnaplus + conckplus + 2. * roots [ j ].
r - conctest - 2. * concsulfate1;
12725 if ( fabs(sumtest) < sumbest ) {
12726 sumbest = fabs(sumtest);
12727 concohminus = conctest;
12742 if ( Istrength < 1.0 ) {
12746 while ( ( fabs(Istrength - Inew) / Istrength ) > 0.10 && counter < 4 ) {
12752 activityCa = ( -Anow *
zCa *
zCa * sqrt(Istrength) ) / ( 1. +
aCa * Bnow * sqrt(Istrength) );
12753 activityCa += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zCa *
zCa * Istrength / sqrt(1000.);
12754 activityCa = exp(activityCa);
12755 activityOH = ( -Anow *
zOH *
zOH * sqrt(Istrength) ) / ( 1. +
aOH * Bnow * sqrt(Istrength) );
12756 activityOH += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zOH *
zOH * Istrength / sqrt(1000.);
12757 activityOH = exp(activityOH);
12758 activityK = ( -Anow *
zK *
zK * sqrt(Istrength) ) / ( 1. +
aK * Bnow * sqrt(Istrength) );
12759 activityK += ( 0.2 - 0.0000417 * Istrength ) * Anow *
zK *
zK * Istrength / sqrt(1000.);
12760 activityK = exp(activityK);
12762 concohminus = conckplus + concnaplus;
12763 if ( (
conccaplus ) > ( 0.1 * ( concohminus ) ) ) {
12767 conccaplus = ( KspCH / ( activityCa * activityOH * activityOH * concohminus * concohminus ) );
12779 if ( syn_old != 2 ) {
12780 test_precip = conckplus * conckplus * activityK * activityK;
12785 printf(
"Syngenite precipitating at cycle %d\n",
icyc);
12787 syngen_change = syn_old = 1;
12789 if ( conckplus > 0.002 ) {
12790 conckplus -= 0.001;
12792 }
else if ( conckplus > 0.0002 ) {
12793 conckplus -= 0.0001;
12807 syngen_change = syn_old = 2;
12810 conckplus += 0.001 *
KperSyn;
12818 }
while ( syngen_change != 0 );
12820 if ( concohminus < ( 0.0000001 ) ) {
12821 concohminus = 0.0000001;
12822 conccaplus = ( KspCH / ( activityCa * activityOH * activityOH * concohminus * concohminus ) );
12825 pH_cur = 14.0 + log10(concohminus * activityOH);
12828 Istrength /= 1000.;
12830 conductivity +=
zOH * concohminus * (
lambdaOH_0 / ( 1. +
GOH * sqrt(Istrength) ) );
12831 conductivity +=
zNa * concnaplus * (
lambdaNa_0 / ( 1. +
GNa * sqrt(Istrength) ) );
12832 conductivity +=
zK * conckplus * (
lambdaK_0 / ( 1. +
GK * sqrt(Istrength) ) );
12838 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);
12844int CemhydMatStatus :: IsSolidPhase(
int phase)
12853void CemhydMatStatus :: burn_phases(
int d1,
int d2,
int d3)
12855 long int icur, inew, ncur, nnew;
12856 int i, j, k, cnt_perc, cnt_tot, kh;
12857 int *nmatx, *nmaty, *nmatz;
12858 int xl, xh, j1, k1, px, py, pz, qx, qy, qz;
12859 int xcn, ycn, zcn, x1, y1, z1, igood;
12860 int *nnewx, *nnewy, *nnewz;
12862 int phase_temp [ 51 ];
12883 for ( k = 0; k <
SYSIZE; k++ ) {
12884 for ( j = 0; j <
SYSIZE; j++ ) {
12885 for ( i = 0; i <
SYSIZE; i++ ) {
12886 newmat [ i ] [ j ] [ k ] =
mic_CSH [ i ] [ j ] [ k ];
12891 ArrPerc [ i ] [ j ] [ k ] = 0;
12897 for ( k = 0; k < 51; k++ ) {
12906 for ( k = 0; k <
SYSIZE; k++ ) {
12907 for ( j = 0; j <
SYSIZE; j++ ) {
12910 for ( kh = 0; kh < 51; kh++ ) {
12911 phase_temp [ kh ] = 0;
12916 while (
last != NULL ) {
12924 px =
cx(i, j, k, d1, d2, d3);
12925 py =
cy(i, j, k, d1, d2, d3);
12926 pz =
cz(i, j, k, d1, d2, d3);
12928 if (
IsSolidPhase(newmat [ px ] [ py ] [ pz ]) == 1 ) {
12930 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
12934 newmat [ px ] [ py ] [ pz ] =
BURNT;
12941 nmatx [ ncur ] = i;
12942 nmaty [ ncur ] = j;
12943 nmatz [ ncur ] = k;
12947 for ( inew = 1; inew <= ncur; inew++ ) {
12948 xcn = nmatx [ inew ];
12949 ycn = nmaty [ inew ];
12950 zcn = nmatz [ inew ];
12952 qx =
cx(xcn, ycn, zcn, d1, d2, d3);
12953 qy =
cy(xcn, ycn, zcn, d1, d2, d3);
12954 qz =
cz(xcn, ycn, zcn, d1, d2, d3);
12957 for ( jnew = 1; jnew <= 6; jnew++ ) {
12988 }
else if ( y1 < 0 ) {
12995 }
else if ( z1 < 0 ) {
13000 if ( ( x1 >= 0 ) && ( x1 <
SYSIZE ) ) {
13001 px =
cx(x1, y1, z1, d1, d2, d3);
13002 py =
cy(x1, y1, z1, d1, d2, d3);
13003 pz =
cz(x1, y1, z1, d1, d2, d3);
13007 if ( (
IsSolidPhase(newmat [ px ] [ py ] [ pz ]) == 1 ) &&
13008 ( newmat [ px ] [ py ] [ pz ] !=
C3S ) &&
13009 ( newmat [ px ] [ py ] [ pz ] !=
C2S ) &&
13010 ( newmat [ px ] [ py ] [ pz ] !=
C3A ) &&
13011 ( newmat [ px ] [ py ] [ pz ] !=
C4AF ) ) {
13015 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
13017 newmat [ px ] [ py ] [ pz ] =
BURNT;
13020 printf(
"error in size of nnew %ld\n", nnew);
13023 nnewx [ nnew ] = x1;
13024 nnewy [ nnew ] = y1;
13025 nnewz [ nnew ] = z1;
13033 ( ( newmat [ px ] [ py ] [ pz ] ==
C3S ) ||
13034 ( newmat [ px ] [ py ] [ pz ] ==
C2S ) ||
13035 ( newmat [ px ] [ py ] [ pz ] ==
C3A ) ||
13036 ( newmat [ px ] [ py ] [ pz ] ==
C4AF ) ) ) {
13037 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
13039 newmat [ px ] [ py ] [ pz ] =
BURNT;
13042 printf(
"error in size of nnew %ld\n", nnew);
13045 nnewx [ nnew ] = x1;
13046 nnewy [ nnew ] = y1;
13047 nnewz [ nnew ] = z1;
13051 else if ( (
micpart [ qx ] [ qy ] [ qz ] ==
micpart [ px ] [ py ] [ pz ] ) &&
13052 ( ( newmat [ px ] [ py ] [ pz ] ==
C3S ) ||
13053 ( newmat [ px ] [ py ] [ pz ] ==
C2S ) ||
13054 ( newmat [ px ] [ py ] [ pz ] ==
C3A ) ||
13055 ( newmat [ px ] [ py ] [ pz ] ==
C4AF ) ) &&
13056 ( (
mic_CSH [ qx ] [ qy ] [ qz ] ==
C3S ) ||
13061 phase_temp [ ( int ) newmat [ px ] [ py ] [ pz ] ]++;
13063 newmat [ px ] [ py ] [ pz ] =
BURNT;
13066 printf(
"error in size of nnew %ld\n", nnew);
13069 nnewx [ nnew ] = x1;
13070 nnewy [ nnew ] = y1;
13071 nnewz [ nnew ] = z1;
13085 for ( icur = 1; icur <= ncur; icur++ ) {
13086 nmatx [ icur ] = nnewx [ icur ];
13087 nmaty [ icur ] = nnewy [ icur ];
13088 nmatz [ icur ] = nnewz [ icur ];
13091 }
while ( nnew > 0 );
13097 for ( j1 = 0; j1 <
SYSIZE; j1++ ) {
13098 for ( k1 = 0; k1 <
SYSIZE; k1++ ) {
13099 px =
cx(xl, j1, k1, d1, d2, d3);
13100 py =
cy(xl, j1, k1, d1, d2, d3);
13101 pz =
cz(xl, j1, k1, d1, d2, d3);
13102 qx =
cx(xh, j1, k1, d1, d2, d3);
13103 qy =
cy(xh, j1, k1, d1, d2, d3);
13104 qz =
cz(xh, j1, k1, d1, d2, d3);
13105 if ( ( newmat [ px ] [ py ] [ pz ] ==
BURNT ) && ( newmat [ qx ] [ qy ] [ qz ] ==
BURNT ) ) {
13109 if ( newmat [ px ] [ py ] [ pz ] ==
BURNT ) {
13110 newmat [ px ] [ py ] [ pz ] =
BURNT + 1;
13113 if ( newmat [ qx ] [ qy ] [ qz ] ==
BURNT ) {
13114 newmat [ qx ] [ qy ] [ qz ] =
BURNT + 1;
13119 if ( igood == 2 ) {
13122 while (
last != NULL ) {
13132 for ( kh = 0; kh < 51; kh++ ) {
13133 phase [ kh ] += phase_temp [ kh ];
13154 for ( kh = 0; kh < 51; kh++ ) {
13155 cnt_perc +=
phase [ kh ];
13165 for ( kh = 1; kh <=
HDCSH; kh++ ) {
13170 cnt_tot +=
count [ kh ];
13174 fprintf(
perc_phases,
"%d %f | ", cnt_tot, (
double ) cnt_perc / cnt_tot);
13175 for ( kh = 0; kh <=
HDCSH; kh++ ) {
13208int CemhydMatStatus :: IsConnected(
int cx,
int cy,
int cz,
int dx,
int dy,
int dz)
13210 int CentPhase, NeighPhase;
13220 if ( ( NeighPhase !=
C3S ) &&
13221 NeighPhase !=
C2S &&
13222 NeighPhase !=
C3A &&
13223 NeighPhase !=
C4AF ) {
13228 else if ( ( CentPhase !=
C3S &&
13229 CentPhase !=
C2S &&
13230 CentPhase !=
C3A &&
13231 CentPhase !=
C4AF ) &&
13232 ( NeighPhase ==
C3S ||
13233 NeighPhase ==
C2S ||
13234 NeighPhase ==
C3A ||
13235 NeighPhase ==
C4AF ) ) {
13243 ( CentPhase ==
C3S ||
13244 CentPhase ==
C2S ||
13245 CentPhase ==
C3A ||
13246 CentPhase ==
C4AF ) &&
13247 ( NeighPhase ==
C3S ||
13248 NeighPhase ==
C2S ||
13249 NeighPhase ==
C3A ||
13250 NeighPhase ==
C4AF ) ) {
13300void CemhydMatStatus :: GenerateConnNumbers()
13554void CemhydMatStatus :: outputImageFilePerc()
13557 char extension [ 10 ];
13560 prefix = (
char * ) malloc(80);
13567 sprintf(extension,
"%04d",
icyc);
13568 strcpy(prefix,
"perc/out5.");
13569 strcat(prefix, extension);
13570 strcat(prefix,
".p.img");
13572 printf(
"Name of percolated output file is %s\n", prefix);
13576 if ( ( perc_img = fopen(prefix,
"w") ) == NULL ) {
13577 printf(
"\nFile %s can not be opened\n", prefix);
13582 for (
int dz = 0; dz <
SYSIZE; dz++ ) {
13583 for (
int dy = 0; dy <
SYSIZE; dy++ ) {
13584 for (
int dx = 0; dx <
SYSIZE; dx++ ) {
13586 fprintf(perc_img,
"%d\n",
ArrPerc [ dx ] [ dy ] [ dz ]);
13593 printf(
"Percolated file %s wrote\n", prefix);
13600void CemhydMatStatus :: WriteUnsortedList(
int px,
int py,
int pz)
13603 if (
last != NULL ) {
13618inline int CemhydMatStatus :: AdjCoord(
int coord)
13624 if ( coord >=
SYSIZE ) {
13632int CemhydMatStatus :: NumSol(
int cx,
int cy,
int cz)
13661 if ( * p_arr ==
CSH ) {
13716void CemhydMatStatus :: nrerror(
const char *error_text)
13718 printf(
"\nNumerical Recipes run-time error...\n");
13719 printf(
"%s\n", error_text);
13720 printf(
"...now exiting to system...\n");
13725float *CemhydMatStatus :: vector(
int nl,
int nh)
13729 v = (
float * ) malloc( (
unsigned ) ( nh - nl + 1 ) *
sizeof(
float ) );
13731 nrerror(
"allocation failure in vector()");
13737int *CemhydMatStatus :: ivector(
int nl,
int nh)
13741 v = (
int * ) malloc( (
unsigned ) ( nh - nl + 1 ) *
sizeof(
int ) );
13743 nrerror(
"allocation failure in ivector()");
13749double *CemhydMatStatus :: dvector(
int nl,
int nh)
13753 v = (
double * ) malloc( (
unsigned ) ( nh - nl + 1 ) *
sizeof(
double ) );
13755 nrerror(
"allocation failure in dvector()");
13764float **CemhydMatStatus :: matrix_cem(
int nrl,
int nrh,
int ncl,
int nch)
13769 m = (
float ** ) malloc( (
unsigned ) ( nrh - nrl + 1 ) *
sizeof(
float * ) );
13771 nrerror(
"allocation failure 1 in matrix_cem()");
13776 for ( i = nrl; i <= nrh; i++ ) {
13777 m [ i ] = (
float * ) malloc( (
unsigned ) ( nch - ncl + 1 ) *
sizeof(
float ) );
13779 nrerror(
"allocation failure 2 in matrix_cem()");
13788double **CemhydMatStatus :: dmatrix(
int nrl,
int nrh,
int ncl,
int nch)
13793 m = (
double ** ) malloc( (
unsigned ) ( nrh - nrl + 1 ) *
sizeof(
double * ) );
13795 nrerror(
"allocation failure 1 in dmatrix()");
13800 for ( i = nrl; i <= nrh; i++ ) {
13801 m [ i ] = (
double * ) malloc( (
unsigned ) ( nch - ncl + 1 ) *
sizeof(
double ) );
13803 nrerror(
"allocation failure 2 in dmatrix()");
13812int **CemhydMatStatus :: imatrix(
int nrl,
int nrh,
int ncl,
int nch)
13816 m = (
int ** ) malloc( (
unsigned ) ( nrh - nrl + 1 ) *
sizeof(
int * ) );
13818 nrerror(
"allocation failure 1 in imatrix()");
13823 for ( i = nrl; i <= nrh; i++ ) {
13824 m [ i ] = (
int * ) malloc( (
unsigned ) ( nch - ncl + 1 ) *
sizeof(
int ) );
13826 nrerror(
"allocation failure 2 in imatrix()");
13837float **CemhydMatStatus :: submatrix(
float **a,
int oldrl,
int oldrh,
int oldcl,
int,
int newrl,
int newcl)
13842 m = (
float ** ) malloc( (
unsigned ) ( oldrh - oldrl + 1 ) *
sizeof(
float * ) );
13844 nrerror(
"allocation failure in submatrix()");
13849 for ( i = oldrl, j = newrl; i <= oldrh; i++, j++ ) {
13850 m [ j ] = a [ i ] + oldcl - newcl;
13868 for ( i = nrh; i >= nrl; i-- ) {
13869 free( (
char * ) ( m [ i ] + ncl ) );
13872 free( (
char * ) ( m + nrl ) );
13879 for ( i = nrh; i >= nrl; i-- ) {
13880 free( (
char * ) ( m [ i ] + ncl ) );
13883 free( (
char * ) ( m + nrl ) );
13890 for ( i = nrh; i >= nrl; i-- ) {
13891 free( (
char * ) ( m [ i ] + ncl ) );
13894 free( (
char * ) ( m + nrl ) );
13899 free( (
char * ) ( b + nrl ) );
13902float **CemhydMatStatus :: convert_matrix(
float *a,
int nrl,
int nrh,
int ncl,
int nch)
13904 int i, j, nrow, ncol;
13907 nrow = nrh - nrl + 1;
13908 ncol = nch - ncl + 1;
13909 m = (
float ** ) malloc( (
unsigned ) ( nrow ) *
sizeof(
float * ) );
13911 nrerror(
"allocation failure in convert_matrix()");
13915 for ( i = 0, j = nrl; i <= nrow - 1; i++, j++ ) {
13916 m [ j ] = a + ncol * i - ncl;
13924 free( (
char * ) ( b + nrl ) );
13946 c.
r = a.
r * b.
r - a.
i * b.
i;
13947 c.
i = a.
i * b.
r + a.
r * b.
i;
13971 if ( fabs(b.
r) >= fabs(b.
i) ) {
13973 den = b.
r + r * b.
i;
13974 c.
r = ( a.
r + r * a.
i ) / den;
13975 c.
i = ( a.
i - r * a.
r ) / den;
13978 den = b.
i + r * b.
r;
13979 c.
r = ( a.
r * r + a.
i ) / den;
13980 c.
i = ( a.
i * r - a.
r ) / den;
13988 float x, y, ans, temp;
13993 }
else if ( y == 0.0 ) {
13995 }
else if ( x > y ) {
13997 ans = x * sqrt(1.0 + temp * temp);
14000 ans = y * sqrt(1.0 + temp * temp);
14010 if ( ( z.
r == 0.0 ) && ( z.
i == 0.0 ) ) {
14019 w = sqrt(x) * sqrt( 0.5 * ( 1.0 + sqrt(1.0 + r * r) ) );
14022 w = sqrt(y) * sqrt( 0.5 * ( r + sqrt(1.0 + r * r) ) );
14025 if ( z.
r >= 0.0 ) {
14027 c.
i = z.
i / ( 2.0 * w );
14029 c.
i = ( z.
i >= 0 ) ? w : -w;
14030 c.
r = z.
i / ( 2.0 * c.
i );
14047void CemhydMatStatus :: CreateHDCSH()
14053void CemhydMatStatus :: PercolateForOutput()
14059double CemhydMatStatus :: GiveWcr()
14068void CemhydMatStatus :: GetInputParams(
char *my_string)
14070 sprintf(my_string,
"%d",
iseed);
14073long CemhydMatStatus :: cx(
int x,
int y,
int z,
int a,
int b,
int c)
14075 return ( 1 - b - c ) * x + ( 1 - a - c ) * y + ( 1 - a - b ) * z;
14079long CemhydMatStatus :: cy(
int x,
int y,
int z,
int a,
int b,
int c)
14081 return ( 1 - a - b ) * x + ( 1 - b - c ) * y + ( 1 - a - c ) * z;
14084long CemhydMatStatus :: cz(
int x,
int y,
int z,
int a,
int b,
int c)
14086 return ( 1 - a - c ) * x + ( 1 - a - b ) * y + ( 1 - b - c ) * z;
14090void CemhydMatStatus :: AnalyticHomogenizationPaste(
double &
E,
double &nu,
int perc_unperc_flag)
14097 double E_CSH_hmg, nu_CSH_hmg;
14103 for ( x = 0; x < 34; x++ ) {
14107 for ( x = 0; x < 51; x++ ) {
14109 if ( x ==
HDCSH ) {
14111 }
else if ( x ==
EMPTYP ) {
14118 if ( perc_unperc_flag == 0 ) {
14129 for ( x = 2; x < 34; x++ ) {
14138 for ( x = 0; x < 34; x++ ) {
14144 PhaseMatrix(0, 1) = 0.001;
14145 PhaseMatrix(0, 2) = 0.499924;
14147 PhaseMatrix(1, 1) = 135.0;
14148 PhaseMatrix(1, 2) = 0.300000;
14150 PhaseMatrix(2, 1) = 130.0;
14151 PhaseMatrix(2, 2) = 0.300000;
14153 PhaseMatrix(3, 1) = 145.0;
14154 PhaseMatrix(3, 2) = 0.300000;
14156 PhaseMatrix(4, 1) = 125.0;
14157 PhaseMatrix(4, 2) = 0.300000;
14159 PhaseMatrix(5, 1) = 30.00;
14160 PhaseMatrix(5, 2) = 0.300000;
14162 PhaseMatrix(6, 1) = 62.90;
14163 PhaseMatrix(6, 2) = 0.300000;
14165 PhaseMatrix(7, 1) = 73.60;
14166 PhaseMatrix(7, 2) = 0.275000;
14168 PhaseMatrix(8, 1) = 72.80;
14169 PhaseMatrix(8, 2) = 0.167000;
14171 PhaseMatrix(9, 1) = 79.60;
14172 PhaseMatrix(9, 2) = 0.309900;
14174 PhaseMatrix(10, 1) = 72.80;
14175 PhaseMatrix(10, 2) = 0.16700;
14177 PhaseMatrix(11, 1) = 72.80;
14178 PhaseMatrix(11, 2) = 0.16700;
14180 PhaseMatrix(12, 1) = 72.80;
14181 PhaseMatrix(12, 2) = 0.16700;
14183 PhaseMatrix(13, 1) = 38.00;
14184 PhaseMatrix(13, 2) = 0.30500;
14185 PhaseMatrix(14, 0) = -10.;
14186 PhaseMatrix(14, 1) = -10.0;
14187 PhaseMatrix(14, 2) = -10.000;
14189 PhaseMatrix(15, 1) = 22.40;
14190 PhaseMatrix(15, 2) = 0.25000;
14192 PhaseMatrix(16, 1) = 22.40;
14193 PhaseMatrix(16, 2) = 0.25000;
14195 PhaseMatrix(17, 1) = 22.40;
14196 PhaseMatrix(17, 2) = 0.25000;
14198 PhaseMatrix(18, 1) = 42.30;
14199 PhaseMatrix(18, 2) = 0.32380;
14201 PhaseMatrix(19, 1) = 22.40;
14202 PhaseMatrix(19, 2) = 0.24590;
14204 PhaseMatrix(20, 1) = 22.40;
14205 PhaseMatrix(20, 2) = 0.24590;
14207 PhaseMatrix(21, 1) = 22.40;
14208 PhaseMatrix(21, 2) = 0.24590;
14210 PhaseMatrix(22, 1) = 42.30;
14211 PhaseMatrix(22, 2) = 0.32380;
14213 PhaseMatrix(23, 1) = 22.40;
14214 PhaseMatrix(23, 2) = 0.24590;
14216 PhaseMatrix(24, 1) = 22.40;
14217 PhaseMatrix(24, 2) = 0.25000;
14219 PhaseMatrix(25, 1) = 30.00;
14220 PhaseMatrix(25, 2) = 0.30000;
14222 PhaseMatrix(26, 1) = 60.00;
14223 PhaseMatrix(26, 2) = 0.30000;
14225 PhaseMatrix(27, 1) = 42.30;
14226 PhaseMatrix(27, 2) = 0.32380;
14228 PhaseMatrix(28, 1) = 79.60;
14229 PhaseMatrix(28, 2) = 0.30990;
14231 PhaseMatrix(29, 1) = 30.00;
14232 PhaseMatrix(29, 2) = 0.30000;
14234 PhaseMatrix(30, 1) = 0.001;
14235 PhaseMatrix(30, 2) = 0.00100;
14246 LevelI(0, 1) = 21.7;
14247 LevelI(0, 2) = 0.24;
14249 LevelI(1, 1) = 29.4;
14250 LevelI(1, 2) = 0.24;
14256 for ( i = 0; i < 31; i++ ) {
14257 PhaseMatrix(i, 0) =
PhaseFrac [ i + 1 ];
14263 PhaseMatrix(14, 1) = CSH_level.
E_hmg;
14264 PhaseMatrix(14, 2) = CSH_level.
nu_hmg;
14266 E_CSH_hmg = CSH_level.
E_hmg;
14267 nu_CSH_hmg = CSH_level.
nu_hmg;
14272 nu = Paste_level.
nu_hmg;
14274 ( void ) E_CSH_hmg;
14275 ( void ) nu_CSH_hmg;
14279void CemhydMatStatus :: 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)
14281 Homogenize Paste_level, Mortar_level, Concrete_level;
14285 FloatMatrix PhaseMatrix(4, 3), Mortar(4, 3), Concrete(4, 3);
14288 PhaseMatrix(0, 1) = E_paste_inp;
14289 PhaseMatrix(0, 2) = nu_paste_inp;
14297 PhaseMatrix(3, 1) = 0.001;
14298 PhaseMatrix(3, 2) = 0.001;
14301 * E_paste = Paste_level.
E_hmg;
14302 * nu_paste = Paste_level.
nu_hmg;
14307 double n_FA, n_CA, vol_ITZ_FA, vol_ITZ_CA;
14312 double vol_tot_mortar = vol_tot_paste +
Vol_FA;
14314 Mortar(0, 0) =
Vol_FA / vol_tot_mortar;
14317 Mortar(1, 0) = vol_ITZ_FA / vol_tot_mortar;
14319 Mortar(1, 2) = * nu_paste;
14320 Mortar(2, 0) = ( vol_tot_paste - vol_ITZ_FA ) / vol_tot_mortar;
14321 Mortar(2, 1) = * E_paste;
14322 Mortar(2, 2) = * nu_paste;
14324 Mortar(3, 1) = 10.;
14325 Mortar(3, 2) = 0.3;
14328 * E_mortar = Mortar_level.
E_hmg;
14329 * nu_mortar = Mortar_level.
nu_hmg;
14332 double vol_tot_concrete = vol_tot_mortar +
Vol_CA;
14334 Concrete(0, 0) =
Vol_CA / vol_tot_concrete;
14337 Concrete(1, 0) = vol_ITZ_CA / vol_tot_concrete;
14339 Concrete(1, 2) = * nu_paste;
14340 Concrete(2, 0) = ( vol_tot_mortar - vol_ITZ_CA ) / vol_tot_concrete;
14341 Concrete(2, 1) = * E_mortar;
14342 Concrete(2, 2) = * nu_mortar;
14343 Concrete(3, 0) = 0.;
14344 Concrete(3, 1) = 10.;
14345 Concrete(3, 2) = 0.3;
14348 E_concrete = Concrete_level.
E_hmg;
14349 nu_concrete = Concrete_level.
nu_hmg;
14353void CemhydMatStatus :: GetInitClinkerPhases(
double &c3s,
double &c2s,
double &c3a,
double &c4af,
double &gypsum,
double &hemi,
double &anh)
14373 TransportMaterialStatus :: updateYourself(tStep);
14382double CemhydMatStatus :: giveAverageTemperature()
14388 if ( !cemhydmat->
eachGP ) {
14399CemhydMatStatus :: printOutputAt(FILE *file,
TimeStep *tStep)
const
14404 TransportMaterialStatus :: printOutputAt(file, tStep);
14405 fprintf(file,
" status {");
14406 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) );
14407 if ( ms->Calculate_elastic_homogenization ) {
14408 fprintf(file,
" EVirginPaste %f NuVirginPaste %f EConcrete %f NuConcrete %f", ms->last_values [ 2 ], ms->last_values [ 3 ], ms->last_values [ 4 ], ms->last_values [ 5 ]);
14412 if ( !cemhydmat->
eachGP ) {
14414 fprintf( file,
" master material %d", cemhydmat->
giveNumber() );
14416 fprintf( file,
" slave of material %d", cemhydmat->
giveNumber() );
14419 fprintf( file,
" independent microstructure %p from material %d",
this, cemhydmat->
giveNumber() );
14422 fprintf(file,
"}\n");
14429CemhydMatStatus :: CemhydMatStatus()
14434 printf(
"Constructor of CemhydMatStatus called\n");
14441CemhydMatStatus :: InitializePy(
const char *inp)
14443 initializeMicrostructure();
14444 readInputFileAndInitialize(inp, (
bool ) 1);
14453 #include <boost/python.hpp>
14454BOOST_PYTHON_MODULE(cemhydmodule)
14457 boost :: python :: class_< oofem :: CemhydMatStatus >(
"cemhydmatstatus")
14458 .def(
"InitializePy", & oofem :: CemhydMatStatus :: InitializePy)
14459 .def(
"MoveCycles", & oofem :: CemhydMatStatus :: MoveCycles)
#define _IFT_CemhydMat_nowarnings
#define _IFT_CemhydMat_scaling
#define _IFT_CemhydMat_eachgp
#define _IFT_CemhydMat_capacitytype
#define _IFT_CemhydMat_densitytype
#define _IFT_CemhydMat_inputFileName
#define _IFT_CemhydMat_conductivitytype
#define _IFT_CemhydMat_reinforcementDegree
#define REGISTER_Material(class)
void QueryNumAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, int &val)
int movech(int xcur, int ycur, int zcur, int finalstep, float nucprob)
double Vol_cement_clinker_gypsum
int extfreidel(int xpres, int ypres, int zpres)
unsigned int * CSH_vicinity
int IsConnected(int cx, int cy, int cz, int dx, int dy, int dz)
int movefh3(int xcur, int ycur, int zcur, int finalstep, float nucprob)
CemhydMatStatus(GaussPoint *gp, CemhydMatStatus *CemStat, CemhydMat *cemhydmat, bool withMicrostructure)
int edgecnt(int xck, int yck, int zck, int ph1, int ph2, int ph3)
void rand3d(int phasein, int phaseout, float xpt)
void makeinert(long int ndesire)
void zroots(fcomplex_cem a[], int m, fcomplex_cem roots[], int polish)
int IsSolidPhase(int phase)
int countem(int xp, int yp, int zp, int phin)
long cy(int x, int y, int z, int a, int b, int c)
void laguer(fcomplex_cem a[], int m, fcomplex_cem *x, float eps, int polish)
double Mass_cement_concrete
void nrerror(const char *error_text)
int burn3d(int npix, int d1, int d2, int d3)
int movecacl2(int xcur, int ycur, int zcur, int finalstep)
void addrand(int randid, long int nneed)
GaussPoint * gp
Stores GP of the CemhydMatStatus.
void hydrate(int fincyc, int stepmax, float chpar1, float chpar2, float hgpar1, float hgpar2, float fhpar1, float fhpar2, float gypar1, float gypar2)
double GivePower(double GiveTemp, double TargTime)
int movec3a(int xcur, int ycur, int zcur, int finalstep, float nucprob)
void alloc_int_3D(int ***(&mask), long SYSIZE)
int chksph(int xin, int yin, int zin, int radd, int wflg, int phasein, int phase2)
int movecaco3(int xcur, int ycur, int zcur, int finalstep)
fcomplex_cem Csub(fcomplex_cem a, fcomplex_cem b)
int loccsh(int xcur, int ycur, int zcur, int extent)
fcomplex_cem Csqrt(fcomplex_cem z)
void drawfloc(int xin, int yin, int zin, int radd, int phasein, int phase2)
int Calculate_elastic_homogenization
Flag to proceed percolation filtering and elastic homogenization.
fcomplex_cem ComplexCemhyd(float re, float im)
int movegyp(int xcur, int ycur, int zcur, int finalstep)
double IPVolume
Volume associated to master IP of one CemhydMat.
int movec4a(int xcur, int ycur, int zcur, int finalstep, float nucprob)
double giveAverageTemperature(void)
void QueryStringAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, char *chars)
void WriteUnsortedList(int px, int py, int pz)
void extc3ah6(int xpres, int ypres, int zpres)
int surfpix(int xin, int yin, int zin)
double GiveDoHActual(void)
Return degree of hydration of the receiver.
void setAverageTemperatureVolume(double temperature, double volume)
Auxiliary function for temperature averaging over GPs.
void outputImageFilePerc(void)
int gsphere(int numgen, long int *numeach, int *sizeeach, int *pheach)
long cz(int x, int y, int z, int a, int b, int c)
void disrealnew(double GiveTemp, double TargTime, int flag)
void alloc_shortint_3D(short int ***(&mic), long SYSIZE)
float Cabs(fcomplex_cem z)
struct percolatedpath * current
fcomplex_cem Cdiv(fcomplex_cem a, fcomplex_cem b)
void dealloc_shortint_3D(short int ***(&mic), long SYSIZE)
double * last_values
Array for storing temporary values (elastic properties etc.).
int icyc
Cycle of celular automata.
void outputImageFileUnperc(char ***m)
void alloc_long_3D(long ***(&mic), long SYSIZE)
void PercolateForOutput(void)
double averageTemperature
Average temperature through integration points.
double init_material_time
Inital material time for growing problems.
double PartHeat
The last incremental heat returned from a GP.
int chckedge(int xck, int yck, int zck)
void extfh3(int xpres, int ypres, int zpres)
void extgyps(int xpres, int ypres, int zpres)
void AnalyticHomogenizationPaste(double &E, double &nu, int perc_unperc_flag)
struct percolatedpath * last
void passone(int low, int high, int cycid, int cshexflag)
int moveanh(int xcur, int ycur, int zcur, int finalstep, float nucprgyp)
int movecas2(int xcur, int ycur, int zcur, int finalstep)
void extafm(int xpres, int ypres, int zpres)
void dealloc_int_3D(int ***(&mask), long SYSIZE)
void disrealnew_init(void)
double Vol_entrained_entrapped_air
void extslagcsh(int xpres, int ypres, int zpres)
double Concrete_thermal_conductivity
int movehem(int xcur, int ycur, int zcur, int finalstep, float nucprgyp)
void dealloc_long_3D(long ***(&mic), long SYSIZE)
void CSHbox(unsigned int *CSH_vicinity)
int extstrat(int xpres, int ypres, int zpres)
int moveettr(int xcur, int ycur, int zcur, int finalstep)
fcomplex_cem Cmul(fcomplex_cem a, fcomplex_cem b)
int countbox(int boxsize, int qx, int qy, int qz)
int movepix(int ntomove, int ph1, int ph2)
int moveone(int *xloc, int *yloc, int *zloc, int *act, int sumold)
void sysinit(int ph1, int ph2)
double MoveCycles(double GiveTemp, int cycles)
void alloc_double_3D(double ***(&mic), long SYSIZE)
void initializeMicrostructure(void)
int movecsh(int xcur, int ycur, int zcur, int finalstep, int cycorig)
int chkfloc(int xin, int yin, int zin, int radd)
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)
long int target_anhydrite
double Concrete_bulk_density
void alloc_char_3D(char ***(&mic), long SYSIZE)
void dealloc_double_3D(double ***(&mic), long SYSIZE)
void dealloc_char_3D(char ***(&mic), long SYSIZE)
int readInputFileAndInitialize(const char *inp, bool generateMicrostructure)
fcomplex_cem Cadd(fcomplex_cem a, fcomplex_cem b)
void burn_phases(int d1, int d2, int d3)
fcomplex_cem RCmul(float x, fcomplex_cem a)
void extpozz(int xpres, int ypres, int zpres)
void sysscan(int ph1, int ph2)
int NumSol(int cx, int cy, int cz)
int countboxc(int boxsize, int qx, int qy, int qz)
int extettr(int xpres, int ypres, int zpres, int etype)
long cx(int x, int y, int z, int a, int b, int c)
int moveas(int xcur, int ycur, int zcur, int finalstep)
FloatArray scaling
Array containing scaling factors for density, conductivity and capacity.
virtual int giveCycleNumber(GaussPoint *gp)
Returns cycle number at the closest cycle after the target time.
int conductivityType
Use different methods to evaluate material parameters.
virtual double giveConcreteCapacity(GaussPoint *gp, TimeStep *tStep) const
Returns concrete thermal capacity depending on chosen type.
virtual double giveConcreteDensity(GaussPoint *gp, TimeStep *tStep) const
Returns concrete density depending on chosen type.
int eachGP
Assign a separate microstructure in each integration point.
virtual void storeWeightTemperatureProductVolume(Element *element, TimeStep *tStep)
Store temperatures multiplied with volume around GPs - need before temperature averaging.
virtual double giveTimeOfCycle(GaussPoint *gp)
Returns time of the CEMHYD3D at the first cycle after the target time.
std::string XMLfileName
XML input file name for CEMHYD3D.
double giveIsotropicConductivity(GaussPoint *gp, TimeStep *tStep) const override
Returns concrete heat conductivity depending on chosen type.
virtual double giveDoHActual(GaussPoint *gp)
Returns DoH of the closest CEMHYD3D cycle after the target time.
virtual void clearWeightTemperatureProductVolume(Element *element)
Clear temperatures multiplied with volume around GPs - need before temperature averaging.
int reinforcementDegree
Degree of reinforcement, if defined, reinforcement effect for conductivity and capacity is accounted ...
virtual void averageTemperature()
Perform averaging on a master CemhydMatStatus.
IntArray nowarnings
Array containing warnings supression for density, conductivity, capacity, high temperature.
CemhydMatStatus * MasterCemhydMatStatus
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
virtual double computeVolumeAround(GaussPoint *gp)
double nu_hmg
Effective Poisson's ratio.
void herveZaoui(FloatMatrix &PhaseMatrix)
void selfConsistent(FloatMatrix &PhaseMatrix)
double E_hmg
Effective Young's modulus or the lower bound.
void moriTanaka(FloatMatrix &PhaseMatrix, int refRow)
IsotropicHeatTransferMaterial(int n, Domain *d)
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
double giveTargetTime()
Returns target time.
double giveField() const
Return last field.
TransportMaterialStatus(GaussPoint *g)
#define OOFEM_WARNING(...)
void free_convert_matrix(float **b, int nrl)
void free_submatrix(float *b, int nrl)
void free_dvector(double *v, int nl)
void free_matrix(float **m, int nrl, int nrh, int ncl)
void free_dmatrix(double **m, int nrl, int nrh, int ncl)
void free_imatrix(int **m, int nrl, int nrh, int ncl)
void free_ivector(int *v, int nl)
double sum(const FloatArray &x)
struct oofem::FCOMPLEX fcomplex_cem
void free_vector(float *v, int nl)
struct cluster * nextpart