OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
cemhydmat.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 /* Majority of this file was developed at the National Institute of
36  * Standards and Technology by employees of the Federal
37  * Government in the course of their official duties. Pursuant
38  * to title 17 Section 105 of the United States Code this
39  * software is not subject to copyright protection and is in
40  * the public domain. CEMHYD3D is an experimental system. NIST
41  * assumes no responsibility whatsoever for its use by other
42  * parties, and makes no guarantees, expressed or implied,
43  * about its quality, reliability, or any other characteristic.
44  * We would appreciate acknowledgement if the software is used.
45  * This software can be redistributed and/or modified freely
46  * provided that any derivative works bear some notice that
47  * they are derived from it, and any modified versions bear
48  * some notice that they have been modified.
49  */
50 
51 #include <cstdio>
52 #include <cstdlib>
53 
54 #include "cemhydmat.h"
55 #include "homogenize.h"
56 #include "mathfem.h"
57 
58 #ifdef __TM_MODULE //OOFEM transport module
59  #include "classfactory.h"
60  #include "domain.h"
61  #include "floatmatrix.h"
62  #include "gausspoint.h"
63 #endif
64 
65 namespace oofem {
66 /* This software was developed at the National Institute of */
67 /* Standards and Technology by employees of the Federal */
68 /* Government in the course of their official duties. Pursuant */
69 /* to title 17 Section 105 of the United States Code this */
70 /* software is not subject to copyright protection and is in */
71 /* the public domain. CEMHYD3D is an experimental system. NIST */
72 /* assumes no responsibility whatsoever for its use by other */
73 /* parties, and makes no guarantees, expressed or implied, */
74 /* about its quality, reliability, or any other characteristic. */
75 /* We would appreciate acknowledgement if the software is used. */
76 /* This software can be redistributed and/or modified freely */
77 /* provided that any derivative works bear some notice that */
78 /* they are derived from it, and any modified versions bear */
79 /* some notice that they have been modified. */
80 /* Modified 3/97 to allow placement of pozzolanic, inert and fly ash particles */
81 /* Modified 9/98 to allow placement of various forms of gypsum */
82 /* Documented version produced 1/00 */
83 /* Modified by smilauer@cml.fsv.cvut.cz to include 1 voxel particles 16.6.2005
84  * Dynamical allocation of memory arrays (possible in input file)
85  */
86 
87 //#define OUTFILES //if defined, output files are generated
88 //#define IMAGEFILES //if defined, output percolated and unpercolated images in each cycle (directories perc/ and unperc/)
89 //#define PRINTF //if defined, printf results simultaneously on screen
90 
91 #ifdef __TM_MODULE //OOFEM transport module
92 
93 REGISTER_Material(CemhydMat);
94 
96 {
97  MasterCemhydMatStatus = NULL;
98 }
99 
101 { }
102 
103 //returns hydration power [W/m3 of concrete]
104 void
106 {
107  double averageTemperature;
108  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
109  val.resize(1);
110 
111  if ( eachGP || ms == MasterCemhydMatStatus ) {
112  averageTemperature = ms->giveAverageTemperature();
113  if ( mode == VM_Total || mode == VM_TotalIntrinsic ) {
114  //for nonlinear solver, return the last value even no time has elapsed
115  if ( tStep->giveTargetTime() != ms->LastCallTime ) {
116  val.at(1) = ms->GivePower( averageTemperature, tStep->giveTargetTime() );
117  } else {
118  val.at(1) = ms->PartHeat;
119  }
120  } else {
121  OOFEM_ERROR( "Undefined mode %s\n", __ValueModeTypeToString(mode) );
122  }
123  } else { //return released heat from the master
124  if ( mode == VM_Total || mode == VM_TotalIntrinsic ) {
126  } else {
127  OOFEM_ERROR( "Undefined mode %s\n", __ValueModeTypeToString(mode) );
128  }
129  }
130 
131  //val.at(1) = 1500;//constant source
132 }
133 
134 
136 {
137  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
138  if ( MasterCemhydMatStatus ) {
140  }
141 
142  return ms->GiveCycNum();
143 }
144 
146 {
147  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
148  if ( MasterCemhydMatStatus ) {
150  }
151 
152  return ms->GiveCycTime();
153 }
154 
155 
156 
158 {
159  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
160  if ( MasterCemhydMatStatus ) {
162  }
163 
164  return ms->GiveDoHActual();
165 }
166 
167 //standard units are [Wm-1K-1]
169 {
170  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
171  double conduct = 0.0;
172 
173  if ( MasterCemhydMatStatus ) {
175  }
176 
177  if ( conductivityType == 0 ) { //given from OOFEM input file
178  conduct = IsotropicHeatTransferMaterial :: give('k', gp, tStep);
179  } else if ( conductivityType == 1 ) { //compute according to Ruiz, Schindler, Rasmussen. Kim, Chang: Concrete temperature modeling and strength prediction using maturity concepts in the FHWA HIPERPAV software, 7th international conference on concrete pavements, Orlando (FL), USA, 2001
180  conduct = IsotropicHeatTransferMaterial :: give('k', gp, tStep) * ( 1.33 - 0.33 * ms->GiveDoHActual() );
181  } else {
182  OOFEM_ERROR("Unknown conductivityType %d\n", conductivityType);
183  }
184 
185 
186  //Parallel Voigt model, 20 W/m/K for steel
187  conduct = conduct * ( 1. - this->reinforcementDegree ) + 20. * this->reinforcementDegree;
188 
189  if ( !this->nowarnings.at(2) && ( conduct < 0.3 || conduct > 5 ) ) {
190  OOFEM_WARNING("Weird concrete thermal conductivity %f W/m/K\n", conduct);
191  }
192 
193  conduct *= this->scaling.at(2);
194 
195  return conduct;
196 }
197 
198 //normally it returns J/kg/K of concrete
200 {
201  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
202  double capacityConcrete = 0.0;
203 
204  if ( MasterCemhydMatStatus ) {
206  }
207 
208  if ( capacityType == 0 ) { //given from OOFEM input file
209  capacityConcrete = IsotropicHeatTransferMaterial :: give('c', gp, tStep);
210  } else if ( capacityType == 1 ) { //compute from CEMHYD3D according to Bentz
211  capacityConcrete = ms->computeConcreteCapacityBentz();
212  } else if ( capacityType == 2 ) { //compute from CEMHYD3D directly
213  capacityConcrete = 1000 * ms->GiveCp();
214  } else {
215  OOFEM_ERROR("Unknown capacityType %d\n", capacityType);
216  }
217 
218  //Parallel Voigt model, 500 J/kg/K for steel
219  capacityConcrete = capacityConcrete * ( 1. - this->reinforcementDegree ) + 500. * this->reinforcementDegree;
220 
221  if ( !this->nowarnings.at(3) && ( capacityConcrete < 500 || capacityConcrete > 2000 ) ) {
222  OOFEM_WARNING("Weird concrete heat capacity %f J/kg/K\n", capacityConcrete);
223  }
224 
225  capacityConcrete *= this->scaling.at(3);
226 
227  return capacityConcrete;
228 }
229 
231 {
232  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
233  double concreteBulkDensity = 0.0;
234 
235  if ( MasterCemhydMatStatus ) {
237  }
238 
239  if ( densityType == 0 ) { //get from OOFEM input file
240  concreteBulkDensity = IsotropicHeatTransferMaterial :: give('d', gp, tStep);
241  } else if ( densityType == 1 ) { //get from XML input file
242  concreteBulkDensity = ms->GiveDensity();
243  } else {
244  OOFEM_ERROR("Unknown densityType %d\n", densityType);
245  }
246 
247  //Parallel Voigt model, 7850 kg/m3 for steel
248  concreteBulkDensity = concreteBulkDensity * ( 1. - this->reinforcementDegree ) + 7850. * this->reinforcementDegree;
249 
250  if ( !this->nowarnings.at(1) && ( concreteBulkDensity < 1000 || concreteBulkDensity > 4000 ) ) {
251  OOFEM_WARNING("Weird concrete density %f kg/m3\n", concreteBulkDensity);
252  }
253 
254  concreteBulkDensity *= this->scaling.at(1);
255 
256  return concreteBulkDensity;
257 }
258 
259 double
261 {
262  if ( mode == Capacity ) {
263  return ( giveConcreteCapacity(gp, tStep) * giveConcreteDensity(gp, tStep) );
264  } else if ( mode == IntSource ) { //for nonlinear solver, return dHeat/dTemperature
265  TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
266  //it suffices to compute derivative of Arrhenius equation
267  //double actualTemperature = status->giveTempField().at(1);
268  double lastEquilibratedTemperature = status->giveField().at(1);
269  //double dt = tStep->giveTimeIncrement();
270  double krate, EaOverR, val;
271  CemhydMatStatus *ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
272  if ( MasterCemhydMatStatus ) {
274  }
275 
276  EaOverR = 1000. * ms->E_act / 8.314;
277 
278  if ( ms->icyc > 1 ) {
279  krate = exp( -EaOverR * ( 1. / ( ms->temp_cur + 273.15 ) - 1. / ( lastEquilibratedTemperature + 273.15 ) ) );
280  //use PartHeat from the last cycle as a corrector tangent, at least one cycle has elapsed
281  // if( fabs(3600*ms->time_cur - ms->PrevHydrTime) > 1.e-3 ){
282  // power = ms->heat_new-ms->PrevCycHeat / (3600*ms->time_cur - ms->PrevHydrTime);//[J/s = W] per gram of cement
283  // power *= 1000 * ms->Mass_cement_concrete;//W/m3/s
284  // } else {
285  // power = 1.e-6;
286  // }
287  } else {
288  krate = 1.;
289  }
290 
291  val = EaOverR * krate / ( ms->temp_cur + 273.15 ) / ( ms->temp_cur + 273.15 );
292 
293  return val;
294  } else {
295  OOFEM_ERROR("unknown mode (%s)\n", __MatResponseModeToString(mode) );
296  }
297 
298  return 0.;
299 }
300 
301 
302 int
304 {
305  CemhydMatStatus *ms;
306  if ( MasterCemhydMatStatus ) {
308  } else {
309  ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
310  }
311 
312  if ( type == IST_HydrationDegree ) {
313  answer.resize(1);
314  answer.at(1) = this->giveDoHActual(gp);
315  return 1;
316  } else if ( type == IST_Density ) {
317  answer.resize(1);
318  answer.at(1) = this->giveConcreteDensity(gp,tStep);
319  return 1;
320  } else if ( type == IST_ThermalConductivityIsotropic ) {
321  answer.resize(1);
322  answer.at(1) = this->giveIsotropicConductivity(gp, tStep);
323  return 1;
324  } else if ( type == IST_HeatCapacity ) {
325  answer.resize(1);
326  answer.at(1) = this->giveConcreteCapacity(gp,tStep);
327  return 1;
328  } else if ( type == IST_AverageTemperature ) {
329  answer.resize(1);
330  answer.at(1) = ms->giveAverageTemperature();
331  return 1;
332  } else if ( type == IST_YoungModulusVirginPaste ) {
333  answer.resize(1);
334  answer.at(1) = ms->last_values [ 2 ];
335  return 1;
336  } else if ( type == IST_PoissonRatioVirginPaste ) {
337  answer.resize(1);
338  answer.at(1) = ms->last_values [ 3 ];
339  return 1;
340  } else if ( type == IST_YoungModulusConcrete ) {
341  answer.resize(1);
342  answer.at(1) = ms->last_values [ 4 ];
343  return 1;
344  } else if ( type == IST_PoissonRatioConcrete ) {
345  answer.resize(1);
346  answer.at(1) = ms->last_values [ 5 ];
347  return 1;
348  } else {
349  return TransportMaterial :: giveIPValue(answer, gp, type, tStep);
350  }
351 }
352 
353 int
355 {
356  for ( GaussPoint *gp: *element->giveDefaultIntegrationRulePtr() ) {
357  CemhydMatStatus *ms;
358  if ( !MasterCemhydMatStatus && !eachGP ) {
359  ms = new CemhydMatStatus(1, domain, gp, NULL, this, 1);
361  } else if ( eachGP ) {
362  ms = new CemhydMatStatus(1, domain, gp, MasterCemhydMatStatus, this, 1);
363  } else {
364  ms = new CemhydMatStatus(1, domain, gp, NULL, this, 0);
365  }
366 
367 // if(!gp->giveMaterialStatus()){
368  gp->setMaterialStatus( ms, this->giveNumber() );
369 // }
370  }
371 
372  return 1;
373 }
374 
376 {
377  CemhydMatStatus *ms;
378 
379  for ( GaussPoint *gp: *element->giveDefaultIntegrationRulePtr() ) {
380  ms = static_cast< CemhydMatStatus * >( this->giveStatus(gp) );
381  ms->setAverageTemperatureVolume(0.0, 0.0);
382  }
383 }
384 
386 {
387  FloatArray vecTemperature;
388 
389  if ( !eachGP ) {
390  for ( GaussPoint *gp: *element->giveDefaultIntegrationRulePtr() ) {
391  //when more GPs are lumped to a master GP
392  double dV = element->computeVolumeAround(gp);
393  element->giveIPValue(vecTemperature, gp, IST_Temperature, tStep);
395  }
396  }
397 }
398 
400 {
401  //printf("%f ", MasterCemhydMatStatus->giveAverageTemperature());
402  if ( !eachGP ) {
404  }
405 }
406 
408 {
409  IRResultType result; // Required by IR_GIVE_FIELD macro
410  castingTime = 0.;
411 
412  result = IsotropicHeatTransferMaterial :: initializeFrom(ir); //read d,k,c
413  if ( result != IRRT_OK ) return result;
414 
415  conductivityType = 0;
416  capacityType = 0;
417  densityType = 0;
418  eachGP = 0;
419  nowarnings.resize(4);
420  nowarnings.zero();
421  scaling.resize(3);
422  for ( int i = 1; i <= scaling.giveSize(); i++ ) {
423  scaling.at(i) = 1.;
424  }
425 
426  reinforcementDegree = 0.;
427  //if you want computation of material properties directly from CEMHYD3D, sum up 1 for density, 2 for conductivity, 4 for capacity
433  if ( nowarnings.giveSize() != 4 ) {
434  OOFEM_ERROR("Incorrect size %d of nowarnings", nowarnings.giveSize() );
435  }
436 
438  if ( scaling.giveSize() != 3 ) {
439  OOFEM_ERROR("Incorrect size %d of scaling", nowarnings.giveSize() );
440  }
441 
444 
445  return IRRT_OK;
446 }
447 
448 
451 {
452  OOFEM_ERROR("Use function CemhydMat :: initMaterial instead");
453  return NULL;
454 }
455 
456 
457 //constructor allowing to copy a microstructure from another CemhydMatStatus
458 //particular instance of CemhydMat in an integration point
459 CemhydMatStatus :: CemhydMatStatus(int n, Domain *d, GaussPoint *gp, CemhydMatStatus *CemStat, CemhydMat *cemhydmat, bool withMicrostructure) : TransportMaterialStatus(n, d, gp)
460 {
461  int i, j, k;
462  PartHeat = 0.;
463  //to be sure, set all pointers to NULL
464  mic = NULL;
465  mic_CSH = NULL;
466  micorig = NULL;
467  micpart = NULL;
468  mask = NULL;
469  ArrPerc = NULL;
470  ConnNumbers = NULL;
471  cshage = NULL;
472  faces = NULL;
473  PhaseFrac = NULL;
474  last_values = NULL;
475  phase = NULL;
476 
477  adiafile = NULL;
478  thfile = NULL;
479  elasfile = NULL;
480  heatfile = NULL;
481  chsfile = NULL;
482  ptmpfile = NULL;
483  movfile = NULL;
484  pHfile = NULL;
485  micfile = NULL;
486  fileperc = NULL;
487  percfile = NULL;
488  disprobfile = NULL;
489  phasfile = NULL;
490  perc_phases = NULL;
491  CSHfile = NULL;
492  infoperc = NULL;
493  infoUnperc = NULL;
494 
495  mic = NULL;
496  micorig = NULL;
497  micpart = NULL;
498  cement = NULL;
499  cemreal = NULL;
500  clust = NULL;
501  mask = NULL;
502  curvature = NULL;
503  mic_CSH = NULL;
504  ArrPerc = NULL;
505  ConnNumbers = NULL;
506  faces = NULL;
507 
508  #ifdef TINYXML
509  xmlFile = NULL;
510  #endif
511 
512  headant = NULL;
513  tailant = NULL;
514  soluble = NULL;
515  creates = NULL;
516  CSH_vicinity = NULL;
517  molarvcsh = NULL;
518  watercsh = NULL;
519  xsph = NULL;
520  ysph = NULL;
521  zsph = NULL;
522  iv = NULL;
523  discount = NULL;
524  count = NULL;
525  disprob = NULL;
526  disbase = NULL;
527  specgrav = NULL;
528  molarv = NULL;
529  heatf = NULL;
530  waterc = NULL;
531  pHeffect = NULL;
532  soluble = NULL;
533  creates = NULL;
534 
535 
536  //set common variables in constructor
537  #ifdef PRINTF
538  printf("Constructor of CemhydMatStatus on GP %p, withMicrostructure %d, copy from CemhydMatStatus %p\n", gp, withMicrostructure, CemStat);
539  fflush(stdout);
540  #endif
541  this->gp = gp;
542  if ( withMicrostructure ) {
543  this->initializeMicrostructure();
544  if ( !CemStat ) {
545  this->readInputFileAndInitialize(cemhydmat->XMLfileName.c_str(), 1);
546  } else { //copy 3D microstructure
547  this->readInputFileAndInitialize(cemhydmat->XMLfileName.c_str(), 0); //read input but do not reconstruct 3D microstructure
548  for ( k = 0; k < SYSIZE; k++ ) {
549  for ( j = 0; j < SYSIZE; j++ ) {
550  for ( i = 0; i < SYSIZE; i++ ) {
551  micpart [ i ] [ j ] [ k ] = CemStat->micpart [ i ] [ j ] [ k ];
552  micorig [ i ] [ j ] [ k ] = CemStat->micorig [ i ] [ j ] [ k ];
553  mic [ i ] [ j ] [ k ] = micorig [ i ] [ j ] [ k ];
554  }
555  }
556  }
557  }
558  }
559 }
560 #endif //__TM_MODULE
561 
562 
564 {
565  icyc = 1; //set the cycle counter
566  time_cur = 0.; //hydration time [h]
567  heat_new = 0.;
568  heat_cf = 0.;
569  Cp_now = 0.900; //initial guess [J/g of concrete], which is later updated
570 
571  averageTemperature = 0.;
572  IPVolume = 0.;
573  init_material_time = 0.;
574 
575  xoff [ 0 ] = 1;
576  xoff [ 1 ] = 0;
577  xoff [ 2 ] = 0;
578  xoff [ 3 ] = -1;
579  xoff [ 4 ] = 0;
580  xoff [ 5 ] = 0;
581  xoff [ 6 ] = 1;
582  xoff [ 7 ] = 1;
583  xoff [ 8 ] = -1;
584  xoff [ 9 ] = -1;
585  xoff [ 10 ] = 0;
586  xoff [ 11 ] = 0;
587  xoff [ 12 ] = 0;
588  xoff [ 13 ] = 0;
589  xoff [ 14 ] = 1;
590  xoff [ 15 ] = 1;
591  xoff [ 16 ] = -1;
592  xoff [ 17 ] = -1;
593  xoff [ 18 ] = 1;
594  xoff [ 19 ] = 1;
595  xoff [ 20 ] = 1;
596  xoff [ 21 ] = 1;
597  xoff [ 22 ] = -1;
598  xoff [ 23 ] = -1;
599  xoff [ 24 ] = -1;
600  xoff [ 25 ] = -1;
601  xoff [ 26 ] = 0;
602 
603  yoff [ 0 ] = 0;
604  yoff [ 1 ] = 1;
605  yoff [ 2 ] = 0;
606  yoff [ 3 ] = 0;
607  yoff [ 4 ] = -1;
608  yoff [ 5 ] = 0;
609  yoff [ 6 ] = 1;
610  yoff [ 7 ] = -1;
611  yoff [ 8 ] = 1;
612  yoff [ 9 ] = -1;
613  yoff [ 10 ] = 1;
614  yoff [ 11 ] = -1;
615  yoff [ 12 ] = 1;
616  yoff [ 13 ] = -1;
617  yoff [ 14 ] = 0;
618  yoff [ 15 ] = 0;
619  yoff [ 16 ] = 0;
620  yoff [ 17 ] = 0;
621  yoff [ 18 ] = 1;
622  yoff [ 19 ] = -1;
623  yoff [ 20 ] = 1;
624  yoff [ 21 ] = -1;
625  yoff [ 22 ] = 1;
626  yoff [ 23 ] = 1;
627  yoff [ 24 ] = -1;
628  yoff [ 25 ] = -1;
629  yoff [ 26 ] = 0;
630 
631  zoff [ 0 ] = 0;
632  zoff [ 1 ] = 0;
633  zoff [ 2 ] = 1;
634  zoff [ 3 ] = 0;
635  zoff [ 4 ] = 0;
636  zoff [ 5 ] = -1;
637  zoff [ 6 ] = 0;
638  zoff [ 7 ] = 0;
639  zoff [ 8 ] = 0;
640  zoff [ 9 ] = 0;
641  zoff [ 10 ] = 1;
642  zoff [ 11 ] = 1;
643  zoff [ 12 ] = -1;
644  zoff [ 13 ] = -1;
645  zoff [ 14 ] = 1;
646  zoff [ 15 ] = -1;
647  zoff [ 16 ] = 1;
648  zoff [ 17 ] = -1;
649  zoff [ 18 ] = 1;
650  zoff [ 19 ] = 1;
651  zoff [ 20 ] = -1;
652  zoff [ 21 ] = -1;
653  zoff [ 22 ] = 1;
654  zoff [ 23 ] = -1;
655  zoff [ 24 ] = 1;
656  zoff [ 25 ] = -1;
657  zoff [ 26 ] = 0;
658 
659  iy = 0; //random generator ran1()
660  LastCallTime = -1.e6; //start from begining (the LastCall in the beginning is -1.e6)
661  alpha_cur = 0.0;
662  alpha_last = 0.; //last degree of hydration
663  LastTargTime = 0.; //stores last call of TargTime
664  TargDoHelas = 0.; //stores DoH at which to perform analytic homogenization
665  //definitions originally from CemhydMat.h
666 
667  //for generation of microstructures
668  CEM = 100; /* and greater */
669  CEMID = 1; /* phase identifier for cement */
670  C2SID = 2; /* phase identified for C2S cement */
671  GYPID = 5; /* phase identifier for gypsum */
672  HEMIHYDRATE = 6; /* phase identifier for hemihydrate */
673  POZZID = 8; /* phase identifier for pozzolanic material - REACTIVE */
674  INERTID = 9; /* phase identifier for inert material - UNREACTIVE */
675  SLAGID = 10; /* phase identifier for slag - REACTIVE */
676  AGG = 28; /* phase identifier for flat aggregate - UNREACTIVE */
677  FLYASH = 30; /* phase identifier for all fly ash components - UNREACTIVE*/
678 
679  //for hydration part
680  POROSITY = 0;
681  C3S = 1;
682  C2S = 2;
683  C3A = 3;
684  C4AF = 4;
685  GYPSUM = 5;
686  HEMIHYD = 6;
687  ANHYDRITE = 7;
688  POZZ = 8;
689  INERT = 9;
690  SLAG = 10;
691  ASG = 11; /* aluminosilicate glass */
692  CAS2 = 12;
693  CH = 13;
694  CSH = 14;
695  C3AH6 = 15;
696  ETTR = 16;
697  ETTRC4AF = 17; /* Iron-rich stable ettringite phase */
698  AFM = 18;
699  FH3 = 19;
700  POZZCSH = 20;
701  SLAGCSH = 21; /* Slag gel-hydration product */
702  CACL2 = 22;
703  FREIDEL = 23; /* Freidel's salt */
704  STRAT = 24; /* stratlingite (C2ASH8) */
705  GYPSUMS = 25; /* Gypsum formed from hemihydrate and anhydrite */
706  CACO3 = 26;
707  AFMC = 27;
708  INERTAGG = 28;
709  ABSGYP = 29;
710  DIFFCSH = 30;
711  DIFFCH = 31;
712  DIFFGYP = 32;
713  DIFFC3A = 33;
714  DIFFC4A = 34;
715  DIFFFH3 = 35;
716  DIFFETTR = 36;
717  DIFFCACO3 = 37;
718  DIFFAS = 38;
719  DIFFANH = 39;
720  DIFFHEM = 40;
721  DIFFCAS2 = 41;
722  DIFFCACL2 = 42;
723  EMPTYP = 45; /*Empty porosity due to self desiccation*/
724  HDCSH = 46;
725  OFFSET = 50; /*Offset for highlighted potentially soluble pixel*/
726 
727  NEIGHBORS = 26; /* number of neighbors to consider (6, 18, or 26) in dissolution */
728  BoxSize = 1; /*int describing vicinity of CSH*/
729  SolidLimit = 27; /*how many solid phase voxels must be in a box (max. <=(2*BoxSize+1)^3)*/
730  MAXTRIES = 150000; /* maximum number of random tries for sphere placement */
731  MAXCYC_SEAL = 30000; /* Maximum number of cycles of sealed hydration (originally MAXCYC in disrealnew.c */
732  NUMSIZES = 100; /* maximum number of different particle sizes */
733  MAXSPH = 10000; /* maximum number of elements in a spherical template */
734 
735  Cp_pozz = 0.75;
736  Cp_CH = 0.75;
737  Cp_h2o = 4.18; /* Cp for free water */
738  Cp_bh2o = 2.2; /* Cp for bound water */
739  WN = 0.23; /* water bound per gram of cement during hydration */
740  WCHSH = 0.06; /* water imbibed per gram of cement during chemical shrinkage (estimate) */
741  CUBEMIN = 3; /* Minimum cube size for checking pore size */
742 
743  DISBIAS = 30.0; /* Dissolution bias- to change all dissolution rates */
744  DISMIN = 0.001; /* Minimum dissolution for C3S dissolution */
745  DISMIN2 = 0.00025; /* Minimum dissolution for C2S dissolution */
746  DISMINSLAG = 0.0001; /* Minimum dissolution for SLAG dissolution */
747  DISMINASG = 0.0005; /* Minimum dissolution for ASG dissolution */
748  DISMINCAS2 = 0.0005; /* Minimum dissolution for CAS2 dissolution */
749  DISMIN_C3A_0 = 0.002; /* Minimum dissolution for C3A dissolution */
750  DISMIN_C4AF_0 = 0.0005; /* Minimum dissolution for C4AF dissolution */
751 
752  C3AH6GROW = 0.01; /* Probability for C3AH6 growth */
753  CHGROW = 1.0; /* Probability for CH growth */
754  CHGROWAGG = 1.0; /* Probability for CH growth on aggregate surface */
755  ETTRGROW = 0.002; /* Probability for ettringite growth */
756  C3AETTR = 0.001; /* Probability for reaction of diffusing C3A with ettringite */
757  C3AGYP = 0.001; /* Probability for diffusing C3A to react with diffusing gypsum */
758  SOLIDC3AGYP = 0.5; /* Probability of solid C3A to react with diffusing sulfate */
759  SOLIDC4AFGYP = 0.1; /* Probability of solid C4AF to react with diffusing sulfate */
760  PPOZZ = 0.05; /* base probability for pozzolanic reaction */
761  PCSH2CSH = 0.002; /* probability for CSH dissolution */
762  A0_CHSOL = 1.325; /* Parameters for variation of CH solubility with */
763  A1_CHSOL = 0.008162; /* temperature (data from Taylor- Cement Chemistry) */
764 
765  BURNT = 70; /* label for a burnt pixel <255 (in char type arrays) */
766  SIZE2D = 49000; /* size of matrices for holding burning locations */
767 
768  SIZESET = 100000;
769  AGRATE = 0.25; /* Probability of gypsum absorption by CSH */
770  VOLFACTOR = 0.00001; /* dm per pixel Note- dm*dm*dm = Liters */
771  MASSFACTOR = 0.0001; /* cm per pixel - specific gravities in g/cm^3 */
772  MMNa = 22.9898;
773  MMK = 39.102;
774  MMNa2O = 61.979;
775  MMK2O = 94.203;
776  BNa = 0.00031; /* From Taylor paper in liters (31 mL/1000/ 100 g) */
777  BK = 0.00020; /* From Taylor paper in liters (20 mL/1000/ 100 g) */
778  BprimeNa = 0.0030; /* From Taylor paper in liters (3 mL/1000/ 1 g POZZ) */
779  BprimeK = 0.0033; /* From Taylor paper in liters (3.3 mL/1000/ 1 g POZZ) */
780  KspCH25C = 0.00000646;
781  KspGypsum = 0.0000263;
782  KspSyngenite = 0.00000010;
783  SpecgravSyngenite = 2.607; /* Source Taylor, H.F.W., Cement Chemistry */
784  KperSyn = 2.0; /* moles of K+ per mole of syngenite */
785  activeA0 = 0.0366; /* A at 295 K (from Ken) */
786  activeB0 = 0.01035; /* B at 295 K (from Ken) */
787  zCa = 2.;
788  zSO4 = 2.;
789  zOH = 1.;
790  zNa = 1.;
791  zK = 1.;
792  aK = 1.33;
793  aCa = 1.;
794  aOH = 3.;
795  aNa = 3.;
796  aSO4 = 4.5; /* Estimate as S ionic radii + O ionic diameter */
797  lambdaOH_0 = 198.0; /* Units: S cm-cm eq.^(-1) */
798  lambdaNa_0 = 50.1;
799  lambdaK_0 = 73.5;
800  lambdaSO4_0 = 39.5;
801  lambdaCa_0 = 29.5; /* Note that CRC has 60./2 for this */
802  GOH = 0.353; /* Units: (eq.^2 mol/L)^(-0.5) */
803  GK = 0.548;
804  GNa = 0.733;
805  GCa = 0.771;
806  GSO4 = 0.877;
807  cm2perL2m = 0.1; /* Conversion from cm2/Liter to 1/m */
808  EPSS = 6.e-8;
809  MAXIT = 100;
810  EPSP = 2.0e-6;
811  MAXM = 100;
812 
813  /*random generator*/
814  IA = 16807;
815  IM = 2147483647;
816  IQ = 127773;
817  IR = 2836;
818  NTAB = 32;
819  EPS = 1.2E-07;
820 
821  PhaseFrac = new double [ 34 ];
822  last_values = new double [ 6 ];
823  last_values [ 2 ] = 0.001; //Young's modulus of virgin paste
824  last_values [ 3 ] = 0.499924; //Poisson's ratio of virgin paste
825  last_values [ 4 ] = 0.001; //Young's modulus of concrete
826  last_values [ 5 ] = 0.499924; //Poisson's ratio of concrete
827  phase = new long int [ 51 ];
828 
829  CSH_vicinity = new unsigned int [ ( 2 * BoxSize + 1 ) * ( 2 * BoxSize + 1 ) * ( 2 * BoxSize + 1 ) + 1 ];
830  molarvcsh = new float [ MAXCYC_SEAL ];
831  watercsh = new float [ MAXCYC_SEAL ];
832  xsph = new int [ MAXSPH ];
833  ysph = new int [ MAXSPH ];
834  zsph = new int [ MAXSPH ];
835 
836  iv = new int [ NTAB ];
837 
838  discount = new long int [ EMPTYP + 1 ];
839  count = new long int [ HDCSH + 1 ];
840  disprob = new float [ HDCSH + 1 ];
841  disbase = new float [ EMPTYP + 1 ];
842  specgrav = new float [ EMPTYP + 1 ];
843  molarv = new float [ EMPTYP + 1 ];
844  heatf = new float [ EMPTYP + 1 ];
845  waterc = new float [ EMPTYP + 1 ];
846  pHeffect = new float [ EMPTYP + 1 ];
847 
848  //set zero to arrays
849  for ( int i = 0; i <= EMPTYP; i++ ) {
850  heatf [ i ] = 0.;
851  waterc [ i ] = 0.;
852  }
853 
854  for ( int i = 0; i <= HDCSH; i++ ) {
855  disprob [ i ] = 0.;
856  }
857 
858  soluble = new int [ EMPTYP + 1 ];
859  creates = new int [ EMPTYP + 1 ];
860 }
861 
862 
863 
864 // destructor
866 {
867 #ifdef OUTFILES
868  if ( fileperc != NULL ) {
869  fclose(fileperc);
870  }
871 
872  if ( disprobfile != NULL ) {
873  fclose(disprobfile);
874  }
875 
876  //fclose(percfile);
877  if ( heatfile != NULL ) {
878  fclose(heatfile);
879  }
880 
881  if ( pHfile != NULL ) {
882  fclose(pHfile);
883  }
884 
885  if ( phasfile != NULL ) {
886  fclose(phasfile);
887  }
888 
889  if ( perc_phases != NULL ) {
890  fclose(perc_phases);
891  }
892 
893  if ( adiafile != NULL ) {
894  fclose(adiafile);
895  }
896 
897  if ( elasfile != NULL ) {
898  fclose(elasfile);
899  }
900 
901  if ( CSHfile != NULL ) {
902  fclose(CSHfile);
903  }
904 
905  if ( infoperc != NULL ) {
906  fclose(infoperc);
907  }
908 
909  if ( infoUnperc != NULL ) {
910  fclose(infoUnperc);
911  }
912 
913 #endif
914  delete [] PhaseFrac;
915  delete [] last_values;
916  delete [] phase;
917 
918 #ifdef CMLFILE
919  delete F; //delete cmlfile
920 #endif
921 #ifdef TINYXML
922  if ( xmlFile != NULL ) {
923  delete xmlFile;
924  }
925 #endif
926 
927  delete [] CSH_vicinity;
928  delete [] molarvcsh;
929  delete [] watercsh;
930 
931  delete [] xsph;
932  delete [] ysph;
933  delete [] zsph;
934  delete [] iv;
935  delete [] discount;
936  delete [] count;
937  delete [] disprob;
938  delete [] disbase;
939  delete [] specgrav;
940  delete [] molarv;
941  delete [] heatf;
942  delete [] waterc;
943  delete [] pHeffect;
944  delete [] soluble;
945  delete [] creates;
946 
947  /* Eliminate the whole list */
948  struct ants *curant;
949 
950  if ( headant != NULL ) { //if hydration did not start sucessfully
951  while ( headant->nextant != NULL ) {
952  curant = headant->nextant;
953  free(headant);
954  headant = curant;
955 #ifdef PRINTF
956  printf("Deallocating headant\n");
957 #endif
958  }
959  }
960 
961  if ( tailant != NULL ) { //if hydration did not start sucessfully
962  free(tailant);
963  }
964 
965 #ifdef PRINTF
966  printf("Deallocating arrays\n");
967  fflush(stdout);
968 #endif
969 
974  dealloc_int_3D(mask, SYSIZE + 1);
979 }
980 
982 {
983  mic = new char ** [ SYSIZE ];
984  for ( int x = 0; x < SYSIZE; x++ ) {
985  mic [ x ] = new char * [ SYSIZE ];
986  for ( int y = 0; y < SYSIZE; y++ ) {
987  mic [ x ] [ y ] = new char [ SYSIZE ];
988  if ( mic [ x ] [ y ] == NULL ) {
989  printf("Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
990  }
991  }
992  }
993 }
994 
996 {
997  if ( mic != NULL ) {
998  for ( int x = 0; x < SYSIZE; x++ ) {
999  for ( int y = 0; y < SYSIZE; y++ ) {
1000  delete [] mic [ x ] [ y ];
1001  }
1002 
1003  delete [] mic [ x ];
1004  }
1005 
1006  delete [] mic;
1007  }
1008 }
1009 
1011 {
1012  mic = new long ** [ SYSIZE ];
1013  for ( int x = 0; x < SYSIZE; x++ ) {
1014  mic [ x ] = new long * [ SYSIZE ];
1015  for ( int y = 0; y < SYSIZE; y++ ) {
1016  mic [ x ] [ y ] = new long [ SYSIZE ];
1017  if ( mic [ x ] [ y ] == NULL ) {
1018  printf("Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1019  }
1020  }
1021  }
1022 }
1023 
1024 
1026 {
1027  if ( mic != NULL ) {
1028  for ( int x = 0; x < SYSIZE; x++ ) {
1029  for ( int y = 0; y < SYSIZE; y++ ) {
1030  delete [] mic [ x ] [ y ];
1031  }
1032 
1033  delete [] mic [ x ];
1034  }
1035 
1036  delete [] mic;
1037  }
1038 }
1039 
1041 {
1042  mic = new int ** [ SYSIZE ];
1043  for ( int x = 0; x < SYSIZE; x++ ) {
1044  mic [ x ] = new int * [ SYSIZE ];
1045  for ( int y = 0; y < SYSIZE; y++ ) {
1046  mic [ x ] [ y ] = new int [ SYSIZE ];
1047  if ( mic [ x ] [ y ] == NULL ) {
1048  printf("Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1049  }
1050  }
1051  }
1052 }
1053 
1054 
1056 {
1057  if ( mic != NULL ) {
1058  for ( int x = 0; x < SYSIZE; x++ ) {
1059  for ( int y = 0; y < SYSIZE; y++ ) {
1060  delete [] mic [ x ] [ y ];
1061  }
1062 
1063  delete [] mic [ x ];
1064  }
1065 
1066  delete [] mic;
1067  }
1068 }
1069 
1070 void CemhydMatStatus :: alloc_shortint_3D(short int ***( &mic ), long SYSIZE)
1071 {
1072  mic = new short int ** [ SYSIZE ];
1073  for ( int x = 0; x < SYSIZE; x++ ) {
1074  mic [ x ] = new short int * [ SYSIZE ];
1075  for ( int y = 0; y < SYSIZE; y++ ) {
1076  mic [ x ] [ y ] = new short int [ SYSIZE ];
1077  if ( mic [ x ] [ y ] == NULL ) {
1078  printf("Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1079  }
1080  }
1081  }
1082 }
1083 
1084 
1086 {
1087  if ( mic != NULL ) {
1088  for ( int x = 0; x < SYSIZE; x++ ) {
1089  for ( int y = 0; y < SYSIZE; y++ ) {
1090  delete [] mic [ x ] [ y ];
1091  }
1092 
1093  delete [] mic [ x ];
1094  }
1095 
1096  delete [] mic;
1097  }
1098 }
1099 
1101 {
1102  mic = new double ** [ SYSIZE ];
1103  for ( int x = 0; x < SYSIZE; x++ ) {
1104  mic [ x ] = new double * [ SYSIZE ];
1105  for ( int y = 0; y < SYSIZE; y++ ) {
1106  mic [ x ] [ y ] = new double [ SYSIZE ];
1107  if ( mic [ x ] [ y ] == NULL ) {
1108  printf("Cannot allocate memory (file %s, line %d)\n", __FILE__, __LINE__);
1109  }
1110  }
1111  }
1112 }
1113 
1114 
1116 {
1117  if ( mic != NULL ) {
1118  for ( int x = 0; x < SYSIZE; x++ ) {
1119  for ( int y = 0; y < SYSIZE; y++ ) {
1120  delete [] mic [ x ] [ y ];
1121  }
1122 
1123  delete [] mic [ x ];
1124  }
1125 
1126  delete [] mic;
1127  }
1128 }
1129 
1130 #ifdef TINYXML
1131 //functions to read int, double and string with error checking
1132 void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, int &val)
1133 {
1134  int success;
1135  char key [ 256 ];
1136  XMLHandle docHandle = XMLHandle(xmlFile);
1137  XMLElement *elemSelected = docHandle.FirstChildElement("cemhyd").FirstChildElement(elementName).ToElement();
1138  if ( elemSelected == NULL ) {
1139  printf("Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1140  exit(0);
1141  }
1142 
1143  sprintf(key, "key%d", position);
1144  success = elemSelected->QueryIntAttribute(key, & val);
1145  if ( success != XML_SUCCESS ) {
1146  printf("Cannot read int value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1147  exit(0);
1148  }
1149 }
1150 
1151 void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, long int &val)
1152 {
1153  int temp;
1154  QueryNumAttributeExt(xmlFile, elementName, position, temp);
1155  val = static_cast< long int >(temp);
1156 }
1157 
1158 void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *xmlFile, const char *elementName, const char *key, int &val)
1159 {
1160  int success;
1161  XMLHandle docHandle = XMLHandle(xmlFile);
1162  XMLElement *elemSelected = docHandle.FirstChildElement("cemhyd").FirstChildElement(elementName).ToElement();
1163  if ( elemSelected == NULL ) {
1164  printf("Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1165  exit(0);
1166  }
1167 
1168  success = elemSelected->QueryIntAttribute(key, & val);
1169  if ( success != XML_SUCCESS ) {
1170  printf("Cannot read int value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1171  exit(0);
1172  }
1173 }
1174 
1175 
1176 void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, double &val)
1177 {
1178  int success;
1179  char key [ 256 ];
1180  XMLHandle docHandle = XMLHandle(xmlFile);
1181  XMLElement *elemSelected = docHandle.FirstChildElement("cemhyd").FirstChildElement(elementName).ToElement();
1182  if ( elemSelected == NULL ) {
1183  printf("Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1184  exit(0);
1185  }
1186 
1187  sprintf(key, "key%d", position);
1188  success = elemSelected->QueryDoubleAttribute(key, & val);
1189  if ( success != XML_SUCCESS ) {
1190  printf("Cannot read double value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1191  exit(0);
1192  }
1193 }
1194 
1195 void CemhydMatStatus :: QueryNumAttributeExt(XMLDocument *xmlFile, const char *elementName, const char *key, double &val)
1196 {
1197  int success;
1198  XMLHandle docHandle = XMLHandle(xmlFile);
1199  XMLElement *elemSelected = docHandle.FirstChildElement("cemhyd").FirstChildElement(elementName).ToElement();
1200  if ( elemSelected == NULL ) {
1201  printf("Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1202  exit(0);
1203  }
1204 
1205  success = elemSelected->QueryDoubleAttribute(key, & val);
1206  if ( success != XML_SUCCESS ) {
1207  printf("Cannot read double value or attribute %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1208  exit(0);
1209  }
1210 }
1211 
1212 void CemhydMatStatus :: QueryStringAttributeExt(XMLDocument *xmlFile, const char *elementName, int position, char *chars)
1213 {
1214  int success;
1215  char key [ 256 ];
1216  XMLHandle docHandle = XMLHandle(xmlFile);
1217  std :: string str1;
1218  XMLElement *elemSelected = docHandle.FirstChildElement("cemhyd").FirstChildElement(elementName).ToElement();
1219  if ( elemSelected == NULL ) {
1220  printf("Cannot find entry %s, terminating, file %s, line %d\n", elementName, __FILE__, __LINE__);
1221  exit(0);
1222  }
1223 
1224  sprintf(key, "key%d", position);
1225  //success = elemSelected->QueryStringAttribute(key, & str1);
1226  // Since ubuntu/debian is still stuck at 2.5.3, lacking QueryStringAttribute.
1227  // Change with above whenever packages are updated.
1228  const char *cstr = elemSelected->Attribute(key);
1229  if ( cstr ) {
1230  success = XML_SUCCESS;
1231  } else {
1232  success = XML_NO_ATTRIBUTE;
1233  }
1234  if ( success != XML_SUCCESS ) {
1235  printf("Cannot read string value or key %s from the entry %s, terminating, file %s, line %d\n", key, elementName, __FILE__, __LINE__);
1236  exit(0);
1237  }
1238  str1 = std :: string(cstr);
1239 
1240  strcpy( chars, str1.c_str() );
1241 }
1242 #endif
1243 
1244 /* read input parameters in file, use XML or cmlfile construction
1245  * allocate necessary arrays (especially those dependent on SYSIZE)
1246  * returns (1) if generation of particles in the Representative Volume Element (RVE) was unsuccessful
1247  */
1248 int CemhydMatStatus :: readInputFileAndInitialize(const char *inp, bool generateMicrostructure)
1249 {
1250  int read_micr;
1251 #ifdef CMLFILE
1252  F = new cmlfile(inp);
1253  // set number of keywords
1254  F->set_labels(54);
1255  F->set_label_string(0, "Rand_seed_num");
1256  F->set_label_string(1, "Input_img_file");
1257  F->set_label_string(2, "Input_id_file");
1258  F->set_label_string(3, "Saturated_sealed");
1259  F->set_label_string(4, "Induction_time");
1260  F->set_label_string(5, "Ea_cement");
1261  F->set_label_string(6, "Ea_pozz");
1262  F->set_label_string(7, "Ea_slag");
1263  F->set_label_string(8, "Beta");
1264  F->set_label_string(9, "Mass_SCM_FA_CA_inert_frac");
1265  F->set_label_string(10, "Mass_cem");
1266  F->set_label_string(11, "Cp_SCM_FA_CA_inert");
1267  F->set_label_string(12, "Cp_cem");
1268  F->set_label_string(13, "Given_microstructure");
1269  F->set_label_string(14, "Output_initial_microstructure");
1270  F->set_label_string(15, "Output_initial_microstructure_img_file");
1271  F->set_label_string(16, "Output_initial_microstructure_id_file");
1272  F->set_label_string(17, "Cycle_freq_perc_pore");
1273  F->set_label_string(18, "Cycle_freq_perc_sol");
1274  F->set_label_string(19, "Total_sodium");
1275  F->set_label_string(20, "Total_potassium");
1276  F->set_label_string(21, "Readily_soluble_sodium");
1277  F->set_label_string(22, "Readily_soluble_potassium");
1278  F->set_label_string(23, "Diffusion_steps_per_cycle");
1279  F->set_label_string(24, "CH_nucleation_probability");
1280  F->set_label_string(25, "CH_scale_factor");
1281  F->set_label_string(26, "Gypsum_nucleation_probability");
1282  F->set_label_string(27, "Gypsum_scale_factor");
1283  F->set_label_string(28, "C3AH6_nucleation_probability");
1284  F->set_label_string(29, "C3AH6_scale_factor");
1285  F->set_label_string(30, "FH3_nucleation_probability");
1286  F->set_label_string(31, "FH3_scale_factor");
1287  F->set_label_string(32, "Microstructure_size");
1288  F->set_label_string(33, "Adiabatic_conditions");
1289  F->set_label_string(34, "Vol_cement_clinker_gypsum");
1290  F->set_label_string(35, "Vol_cement_SCM");
1291  F->set_label_string(36, "Vol_water");
1292  F->set_label_string(37, "Vol_FA");
1293  F->set_label_string(38, "Vol_CA");
1294  F->set_label_string(39, "Vol_inert_filler");
1295  F->set_label_string(40, "Vol_entrained_entrapped_air");
1296  F->set_label_string(41, "Grain_average_FA");
1297  F->set_label_string(42, "Grain_average_CA");
1298  F->set_label_string(43, "ITZ_thickness");
1299  F->set_label_string(44, "ITZ_Young_red");
1300  F->set_label_string(45, "Young_SCM");
1301  F->set_label_string(46, "Poisson_SCM");
1302  F->set_label_string(47, "Young_FA");
1303  F->set_label_string(48, "Poisson_FA");
1304  F->set_label_string(49, "Young_CA");
1305  F->set_label_string(50, "Poisson_CA");
1306  F->set_label_string(51, "Young_inert");
1307  F->set_label_string(52, "Poisson_inert");
1308  F->set_label_string(53, "Calculate_elastic_homogenization");
1309 
1310  // these keywords with #id will be required
1311  F->require(0);
1312  //F->require( 1 ) ;
1313  //F->require( 2 ) ;
1314  F->require(3);
1315  F->require(4);
1316  F->require(5);
1317  F->require(6);
1318  F->require(7);
1319  F->require(8);
1320  F->require(9);
1321  F->require(10);
1322  F->require(11);
1323  F->require(12);
1324  F->require(13);
1325  F->require(14);
1326  F->require(17);
1327  F->require(18);
1328  F->require(19);
1329  F->require(20);
1330  F->require(21);
1331  F->require(22);
1332  F->require(23);
1333  F->require(24);
1334  F->require(27);
1335  F->require(28);
1336  F->require(29);
1337  F->require(30);
1338  F->require(31);
1339  F->require(32);
1340  F->require(33);
1341  F->require(34);
1342  F->require(35);
1343  F->require(36);
1344  F->require(37);
1345  F->require(38);
1346  F->require(39);
1347  F->require(40);
1348  F->require(41);
1349  F->require(42);
1350  F->require(43);
1351  F->require(44);
1352  F->require(45);
1353  F->require(46);
1354  F->require(47);
1355  F->require(48);
1356  F->require(49);
1357  F->require(50);
1358  F->require(51);
1359  F->require(52);
1360  F->require(53);
1361  // set number and names of sections
1362  F->set_sections(1);
1363  F->set_section_string(0, "CEMHYD_generate_particles");
1364 
1365  F->check_requirements();
1366 
1367  if ( F->error_in_requirements() ) {
1368  printf("Cemhyd input file %s is not complete (file %s, line %d)\n", inp, __FILE__, __LINE__);
1369  exit(0);
1370  }
1371 
1372  F->get_value(0, ( long & )iseed);
1373 #endif
1374 #ifdef TINYXML
1375  xmlFile = new XMLDocument();
1376  countKey = 0;
1377  int errorId = xmlFile->LoadFile(inp);
1378  if ( errorId != XML_SUCCESS ) {
1379  printf("\nError reading XML file %s or nonletter symbols used, error id = %d\n", inp, errorId);
1380  exit(0);
1381  }
1382 
1383  QueryNumAttributeExt(xmlFile, "Rand_seed_num", 0, iseed);
1384  QueryNumAttributeExt(xmlFile, "Microstructure_size", 0, SYSSIZE);
1385  QueryNumAttributeExt(xmlFile, "Given_microstructure", 0, read_micr);
1386 #endif
1387 
1388  nseed = iseed;
1389  seed = ( & nseed );
1390  //printf("iseed = %d", iseed);
1391 
1392  //read SYSSIZE of microstructure and allocate arrays
1393 #ifdef CMLFILE
1394  F->get_value(32, ( long & )SYSSIZE);
1395 #endif
1396  if ( SYSSIZE < 10 ) {
1397  printf("Can not run small microstructure %d (< 10 voxels a side), file %s, line %d)\n", SYSSIZE, __FILE__, __LINE__);
1398  exit(0);
1399  }
1400 
1401  SYSIZE = SYSSIZE;
1402  SYSIZEM1 = ( SYSIZE - 1 ); /* System size -1 */
1403  SYSIZE_POW3 = ( SYSIZE * SYSIZE * SYSIZE );
1404  NPARTC = ( long ) ( 700000 * ( double ) SYSIZE_POW3 / 1000000. );
1405  BURNTG = ( NPARTC > 100 ? NPARTC : 100 );
1406  CUBEMAX = ( SYSIZE > 6 ? 7 : 3 );
1407  DETTRMAX = ( 1200. * SYSIZE_POW3 / 1000000. ); /* Maximum allowed # of ettringite diffusing species */
1408  DGYPMAX = ( 2000. * SYSIZE_POW3 / 1000000. ); /* Maximum allowed # of gypsum diffusing species */
1409  DCACO3MAX = ( 1000 * SYSIZE_POW3 / 1000000. ); /* Maximum allowed # of CaCO3 diffusing species */
1410  DCACL2MAX = ( 2000 * SYSIZE_POW3 / 1000000. ); /* Maximum allowed # of CaCl2 diffusing species */
1411  DCAS2MAX = ( 2000 * SYSIZE_POW3 / 1000000. ); /* Maximum allowed # of CAS2 diffusing species */
1412  CHCRIT = ( 50. * SYSIZE_POW3 / 1000000. ); /* Scale parameter to adjust CH dissolution probability */
1413  C3AH6CRIT = ( 10. * SYSIZE_POW3 / 1000000. ); /* Scale par. to adjust C3AH6 dissolution prob. */
1414  CSHSCALE = ( 70000. * SYSIZE_POW3 / 1000000. ); /*scale factor for CSH controlling induction */
1415  C3AH6_SCALE = ( 2000. * SYSIZE_POW3 / 1000000. ); /*scale factor for C3AH6 controlling induction of aluminates */
1416 
1421  alloc_int_3D(mask, SYSIZE + 1);
1426 
1427 #ifdef CMLFILE
1428  F->get_value(13, ( long & )read_micr);
1429 #endif
1430 
1431  if ( !read_micr && generateMicrostructure == 1 ) { //generate new microstructure
1432  if ( genpartnew() == 1 ) { //read input file for RVE generation, if unsuccessful microstructure generation, return
1433  return 1;
1434  }
1435 
1436 #ifdef PRINTF
1437  printf("MONOPHASE microstructure created\n");
1438 #endif
1439  distrib3d(); //read autocorrelation functions and distribute clinker phases in RVE
1440  }
1441 
1442  readhydrparam(); //read hydration parameters
1443 
1444  return 0;
1445 }
1446 
1447 
1448 /* Random number generator ran1 from Computers in Physics */
1449 /* Volume 6 No. 5, 1992, 522-524, Press and Teukolsky */
1450 /* To generate real random numbers 0.0-1.0 */
1451 /* Should be seeded with a negative integer */
1452 double CemhydMatStatus :: ran1(int *idum)
1453 /* Calls: no routines */
1454 /* Called by: gsphere,makefloc */
1455 {
1456  int j, k;
1457  NDIV = 1.0 / ( 1.0 + ( IM - 1.0 ) / NTAB );
1458  RNMX = ( 1.0 - EPS );
1459  AM = ( 1.0 / IM );
1460 
1461  if ( ( * idum <= 0 ) || ( iy == 0 ) ) {
1462  * idum = ( -* idum > * idum ) ? -* idum : * idum; //MAX(-*idum>*idum);
1463  for ( j = NTAB + 7; j >= 0; j-- ) {
1464  k = * idum / IQ;
1465  * idum = IA * ( * idum - k * IQ ) - IR * k;
1466  if ( * idum < 0 ) {
1467  * idum += IM;
1468  }
1469 
1470  if ( j < NTAB ) {
1471  iv [ j ] = * idum;
1472  }
1473  }
1474 
1475  iy = iv [ 0 ];
1476  }
1477 
1478  k = * idum / IQ;
1479  * idum = IA * ( * idum - k * IQ ) - IR * k;
1480  if ( * idum < 0 ) {
1481  * idum += IM;
1482  }
1483 
1484  j = ( int ) ( iy * NDIV );
1485  iy = iv [ j ];
1486  iv [ j ] = * idum;
1487  return AM * iy < RNMX ? AM * iy : RNMX; //MIN(AM*iy,RNMX);
1488 }
1489 
1490 
1491 /* routine to add a flat plate aggregate in the microstructure */
1493 /* Calls: no other routines */
1494 /* Called by: main program */
1495 {
1496  int ix, iy, iz;
1497  int agglo, agghi;
1498 
1499  /* Be sure aggregate size is an even integer */
1500  do {
1501  printf("Enter thickness of aggregate to place (an even integer) \n");
1502 #ifdef CMLFILE
1503  F->get_next_line_in_section(0, ( long & )aggsize);
1504 #endif
1505 #ifdef TINYXML
1506  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, aggsize);
1507 #endif
1508 
1509  //fscanf(in, "%d",&aggsize);
1510 #ifdef PRINTF
1511  printf("%d\n", aggsize);
1512 #endif
1513  } while ( ( ( aggsize % 2 ) != 0 ) || ( aggsize > ( SYSSIZE - 2 ) ) );
1514 
1515  if ( aggsize != 0 ) {
1516  agglo = ( SYSSIZE / 2 ) - ( ( aggsize - 2 ) / 2 );
1517  agghi = ( SYSSIZE / 2 ) + ( aggsize / 2 );
1518 
1519  /* Aggregate is placed in yz plane */
1520  for ( ix = agglo; ix <= agghi; ix++ ) {
1521  for ( iy = 1; iy <= SYSSIZE; iy++ ) {
1522  for ( iz = 1; iz <= SYSSIZE; iz++ ) {
1523  /* Mark aggregate in both particle and microstructure images */
1524  cement [ ix ] [ iy ] [ iz ] = AGG;
1525  cemreal [ ix ] [ iy ] [ iz ] = AGG;
1526  }
1527  }
1528  }
1529  }
1530 }
1531 
1532 
1533 
1534 /* routine to check or perform placement of sphere of ID phasein */
1535 /* centered at location (xin,yin,zin) of radius radd */
1536 /* wflg=1 check for fit of sphere */
1537 /* wflg=2 place the sphere */
1538 /* phasein and phase2 are phases to assign to cement and cemreal images resp. */
1539 int CemhydMatStatus :: chksph(int xin, int yin, int zin, int radd, int wflg, int phasein, int phase2)
1540 /* Calls: no other routines */
1541 /* Called by: gsphere */
1542 {
1543  int nofits, xp, yp, zp, i, j, k;
1544  float dist, xdist, ydist, zdist, ftmp;
1545 
1546  nofits = 0; /* Flag indicating if placement is possible */
1547  /* Check all pixels within the digitized sphere volume */
1548  for ( i = xin - radd; ( ( i <= xin + radd ) && ( nofits == 0 ) ); i++ ) {
1549  xp = i;
1550  /* use periodic boundary conditions for sphere placement */
1551  if ( xp < 1 ) {
1552  xp += SYSSIZE;
1553  } else if ( xp > SYSSIZE ) {
1554  xp -= SYSSIZE;
1555  }
1556 
1557  ftmp = ( float ) ( i - xin );
1558  xdist = ftmp * ftmp;
1559  for ( j = yin - radd; ( ( j <= yin + radd ) && ( nofits == 0 ) ); j++ ) {
1560  yp = j;
1561  /* use periodic boundary conditions for sphere placement */
1562  if ( yp < 1 ) {
1563  yp += SYSSIZE;
1564  } else if ( yp > SYSSIZE ) {
1565  yp -= SYSSIZE;
1566  }
1567 
1568  ftmp = ( float ) ( j - yin );
1569  ydist = ftmp * ftmp;
1570  for ( k = zin - radd; ( ( k <= zin + radd ) && ( nofits == 0 ) ); k++ ) {
1571  zp = k;
1572  /* use periodic boundary conditions for sphere placement */
1573  if ( zp < 1 ) {
1574  zp += SYSSIZE;
1575  } else if ( zp > SYSSIZE ) {
1576  zp -= SYSSIZE;
1577  }
1578 
1579  ftmp = ( float ) ( k - zin );
1580  zdist = ftmp * ftmp;
1581 
1582  /* Compute distance from center of sphere to this pixel */
1583  dist = sqrt(xdist + ydist + zdist);
1584  if ( ( dist - 0.5 ) <= ( float ) radd ) {
1585  /* Perform placement */
1586  if ( wflg == 2 ) {
1587  cement [ xp ] [ yp ] [ zp ] = phasein;
1588  cemreal [ xp ] [ yp ] [ zp ] = phase2;
1589  }
1590  /* or check placement */
1591  else if ( ( wflg == 1 ) && ( cement [ xp ] [ yp ] [ zp ] != POROSITY ) ) {
1592  nofits = 1;
1593  }
1594  }
1595 
1596  /* Check for overlap with aggregate */
1597  if ( ( wflg == 1 ) && ( ( fabs( xp - ( ( float ) ( SYSSIZE + 1 ) / 2.0 ) ) ) < ( ( float ) aggsize / 2.0 ) ) ) {
1598  nofits = 1;
1599  }
1600  }
1601  }
1602  }
1603 
1604  /* return flag indicating if sphere will fit */
1605  return ( nofits );
1606 }
1607 
1608 
1609 
1610 /* routine to place spheres of various sizes and phases at random */
1611 /* locations in 3-D microstructure */
1612 /* numgen is number of different size spheres to place */
1613 /* numeach holds the number of each size class */
1614 /* sizeeach holds the radius of each size class */
1615 /* pheach holds the phase of each size class */
1616 int CemhydMatStatus :: gsphere(int numgen, long int *numeach, int *sizeeach, int *pheach)
1617 {
1618  /* Calls: makesph, ran1 */
1619  /* Called by: create */
1620  int count, x, y, z, radius, ig, tries, phnow;
1621  long int jg, i;
1622  float testgyp, typegyp;
1623 
1624  /* Generate spheres of each size class in turn (largest first) */
1625  for ( ig = 0; ig < numgen; ig++ ) {
1626  phnow = pheach [ ig ]; /* phase for this class */
1627  radius = sizeeach [ ig ]; /* radius for this class */
1628  /* loop for each sphere in this size class */
1629  for ( jg = 1; jg <= numeach [ ig ]; jg++ ) {
1630  tries = 0;
1631  /* Stop after MAXTRIES random tries */
1632  do {
1633  tries += 1;
1634  /* generate a random center location for the sphere */
1635  x = ( int ) ( ( float ) SYSSIZE * ran1(seed) ) + 1;
1636  y = ( int ) ( ( float ) SYSSIZE * ran1(seed) ) + 1;
1637  z = ( int ) ( ( float ) SYSSIZE * ran1(seed) ) + 1;
1638  /* See if the sphere will fit at x,y,z */
1639  /* Include dispersion distance when checking */
1640  /* to insure requested separation between spheres */
1641  count = chksph(x, y, z, radius + dispdist, 1, npart + CEM, 0);
1642  if ( ( tries > MAXTRIES ) && ( dispdist == 2 ) ) {
1643  tries = 0;
1644  dispdist += 1;
1645  }
1646 
1647  if ( tries > MAXTRIES ) {
1648  printf("Could not place sphere %d after %ld random attempts \n", npart, MAXTRIES);
1649  printf("Skipping this microstructure parameters\n");
1650  for ( i = 1; i <= npart; i++ ) {
1651  free(clust [ i ]);
1652  }
1653 
1654  return ( 1 );
1655  }
1656  } while ( count != 0 );
1657 
1658  /* place the sphere at x,y,z */
1659  npart += 1;
1660  if ( npart >= NPARTC ) {
1661  printf("Too many spheres being generated \n");
1662  printf("User needs to increase value of NPARTC at top of C-code\n");
1663  printf("Skipping this microstructure parameters\n");
1664  return ( 1 );
1665  }
1666 
1667  /* Allocate space for new particle info */
1668  clust [ npart ] = ( struct cluster * ) malloc( sizeof( struct cluster ) );
1669  clust [ npart ]->partid = npart;
1670  clust [ npart ]->clustid = npart;
1671  /* Default to cement placement */
1672  clust [ npart ]->partphase = CEMID;
1673  clust [ npart ]->x = x;
1674  clust [ npart ]->y = y;
1675  clust [ npart ]->z = z;
1676  clust [ npart ]->r = radius;
1677  clusleft += 1;
1678  if ( phnow == 1 ) {
1679  testgyp = ran1(seed);
1680  /* Do not use dispersion distance when placing particle */
1681  if ( ( ( testgyp > probgyp ) && ( ( target_sulfate - n_sulfate ) < ( target_total - n_total ) ) ) || ( n_sulfate > target_sulfate ) || ( volpart [ radius ] > ( target_sulfate - n_sulfate ) ) || ( numeach [ ig ] <= 2 ) ) {
1682  count = chksph(x, y, z, radius, 2, npart + CEM - 1, CEMID);
1683  n_total += volpart [ radius ];
1684  } else {
1685  /* Place particle as gypsum */
1686  typegyp = ran1(seed);
1687  n_total += volpart [ radius ];
1688  n_sulfate += volpart [ radius ];
1689  if ( ( probanh >= 1.0 ) || ( ( typegyp < probanh ) && ( n_anhydrite < target_anhydrite ) && ( volpart [ radius ] <= ( target_anhydrite - n_anhydrite ) ) ) ) {
1690  /* Place particle as anhydrite */
1691  n_anhydrite += volpart [ radius ];
1692  count = chksph(x, y, z, radius, 2, npart + CEM - 1, ANHYDRITE);
1693  clust [ npart ]->partphase = ANHYDRITE;
1694  } else if ( ( ( probanh + probhem ) >= 1.0 ) || ( ( typegyp < ( probanh + probhem ) ) && ( n_hemi < target_hemi ) && ( volpart [ radius ] <= ( target_hemi - n_hemi ) ) ) ) {
1695  /* Place particle as hemihydrate */
1696  n_hemi += volpart [ radius ];
1697  count = chksph(x, y, z, radius, 2, npart + CEM - 1, HEMIHYDRATE);
1699  } else {
1700  count = chksph(x, y, z, radius, 2, npart + CEM - 1, GYPID);
1701  /* Correct phase ID of particle */
1702  clust [ npart ]->partphase = GYPID;
1703  }
1704  }
1705  }
1706  /* place as inert, CaCO3, C2S, slag, or pozzolanic material */
1707  else {
1708  count = chksph(x, y, z, radius, 2, npart + CEM - 1, phnow);
1709  /* Correct phase ID of particle */
1710  clust [ npart ]->partphase = phnow;
1711  }
1712 
1713  clust [ npart ]->nextpart = NULL;
1714  }
1715  }
1716 
1717  //deallocate
1718  for ( i = 1; i <= npart; i++ ) {
1719  free(clust [ i ]);
1720  //printf("Dealloc clust %ld of %d", i,npart);
1721  }
1722 
1723  return ( 0 );
1724 }
1725 
1726 
1727 /* routine to obtain user input and create a starting microstructure */
1729 {
1730  /* Calls: gsphere */
1731  /* Called by: main program */
1732  int numsize;
1733  int *sphrad, *sphase;
1734  long int *sphnum;
1735  long int inval1;
1736  int isph, inval;
1737 
1738  sphrad = new int [ NUMSIZES ];
1739  sphase = new int [ NUMSIZES ];
1740  sphnum = new long int [ NUMSIZES ];
1741 
1742 
1743  do {
1744 #ifdef PRINTF
1745  printf("Enter number of different size spheres to use(max. is %d) \n", NUMSIZES);
1746 #endif
1747 #ifdef CMLFILE
1748  F->get_next_line_in_section(0, ( long & )numsize);
1749 #endif
1750 #ifdef TINYXML
1751  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, numsize);
1752 #endif
1753  //fscanf(in, "%d",&numsize);
1754 #ifdef PRINTF
1755  printf("%d \n", numsize);
1756 #endif
1757  } while ( ( numsize > NUMSIZES ) || ( numsize < 0 ) );
1758 
1759  do {
1760 #ifdef PRINTF
1761  printf("Enter dispersion factor (separation distance in pixels) for spheres (0-2) \n");
1762  printf("0 corresponds to totally random placement \n");
1763 #endif
1764 #ifdef CMLFILE
1765  F->get_next_line_in_section(0, ( long & )dispdist);
1766 #endif
1767 #ifdef TINYXML
1768  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, dispdist);
1769 #endif
1770  //fscanf(in, "%d",&dispdist);
1771 #ifdef PRINTF
1772  printf("%d \n", dispdist);
1773 #endif
1774  } while ( ( dispdist < 0 ) || ( dispdist > 2 ) );
1775 
1776  do {
1777 #ifdef PRINTF
1778  printf("Enter probability for gypsum particles on a random particle basis (0.0-1.0) \n");
1779 #endif
1780 #ifdef CMLFILE
1781  F->get_next_line_in_section(0, probgyp);
1782 #endif
1783 #ifdef TINYXML
1784  QueryNumAttributeExt(xmlFile, "Generate_microstructure", "dihydrate", probgyp);
1785 #endif
1786  //fscanf(in, "%f",&probgyp);
1787 #ifdef PRINTF
1788  printf("%f \n", probgyp);
1789 #endif
1790  } while ( ( probgyp < 0.0 ) || ( probgyp > 1.0 ) );
1791 
1792  do {
1793 #ifdef PRINTF
1794  printf("Enter probabilities for hemihydrate and anhydrite forms of gypsum (0.0-1.0) \n");
1795 #endif
1796 #ifdef CMLFILE
1797  F->get_next_line_in_section(0, probhem);
1798  F->get_next_line_in_section(0, probanh);
1799 #endif
1800 #ifdef TINYXML
1801  QueryNumAttributeExt(xmlFile, "Generate_microstructure", "hemihydrate", probhem);
1802  QueryNumAttributeExt(xmlFile, "Generate_microstructure", "anhydrite", probanh);
1803 #endif
1804  //fscanf(in, "%f %f",&probhem,&probanh);
1805 #ifdef PRINTF
1806  printf("%f %f\n", probhem, probanh);
1807 #endif
1808  } while ( ( probhem < 0.0 ) || ( probhem > 1.0 ) || ( probanh < 0.0 ) || ( probanh > 1.0 ) || ( ( probanh + probhem ) > 1.001 ) );
1809 
1810  if ( ( numsize > 0 ) && ( numsize < ( NUMSIZES + 1 ) ) ) {
1811 #ifdef PRINTF
1812  printf("Enter number, radius, and phase ID for each sphere class (largest radius 1st) \n");
1813  printf("Phases are %d- Cement and (random) calcium sulfate, %d- C2S, %d- Gypsum, %d- hemihydrate %d- anhydrite %d- Pozzolanic, %d- Inert, %d- Slag, %d- CaCO3 %d- Fly Ash \n", CEMID, C2SID, GYPID, HEMIHYDRATE, ANHYDRITE, POZZID, INERTID, SLAGID, CACO3, FLYASH);
1814 #endif
1815  /* Obtain input for each size class of spheres */
1816  for ( isph = 0; isph < numsize; isph++ ) {
1817 #ifdef PRINTF
1818  printf("Enter number of spheres of class %d \n", isph + 1);
1819 #endif
1820 #ifdef CMLFILE
1821  F->get_next_line_in_section(0, inval1);
1822 #endif
1823 #ifdef TINYXML
1824  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, inval1);
1825  //inval1 = static_cast<long int>(inval);
1826 #endif
1827  //fscanf(in, "%ld",&inval1);
1828 #ifdef PRINTF
1829  printf("%ld \n", inval1);
1830 #endif
1831  sphnum [ isph ] = inval1;
1832 
1833  // do{
1834 #ifdef PRINTF
1835  printf("Enter radius of spheres of class %d \n", isph + 1);
1836  printf("(Integer <=%d please) \n", SYSSIZE / 3);
1837 #endif
1838 #ifdef CMLFILE
1839  F->get_next_line_in_section(0, ( long & )inval);
1840 #endif
1841 #ifdef TINYXML
1842  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, inval);
1843 #endif
1844  //fscanf(in, "%d",&inval);
1845  if ( inval > ( SYSSIZE / 3 ) ) {
1846  printf("Given radius %d exceeded maximum radius of %d, terminating\n", inval, SYSSIZE / 3);
1847  exit(0);
1848  }
1849 
1850 #ifdef PRINTF
1851  printf("%d \n", inval);
1852 #endif
1853  //} while ((inval<0)||(inval>(SYSSIZE/3)));
1854 
1855  sphrad [ isph ] = inval;
1856  do {
1857 #ifdef PRINTF
1858  printf("Enter phase of spheres of class %d \n", isph + 1);
1859 #endif
1860 #ifdef CMLFILE
1861  F->get_next_line_in_section(0, ( long & )inval);
1862 #endif
1863 #ifdef TINYXML
1864  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, inval);
1865 #endif
1866  //fscanf(in, "%d",&inval);
1867 #ifdef PRINTF
1868  printf("%d \n", inval);
1869 #endif
1870  } while ( ( inval != CEMID ) && ( inval != C2SID ) && ( inval != GYPID ) && ( inval != HEMIHYDRATE ) && ( inval != ANHYDRITE ) && ( inval != POZZID ) && ( inval != INERTID ) && ( inval != SLAGID ) && ( inval != FLYASH ) && ( inval != CACO3 ) && ( inval != AGG ) && ( inval != ASG ) );
1871 
1872  sphase [ isph ] = inval;
1873  if ( inval == CEMID ) {
1874  target_total += sphnum [ isph ] * volpart [ sphrad [ isph ] ];
1875  }
1876  }
1877 
1878  /* Determine target pixel counts for calcium sulfate forms */
1879  target_sulfate = ( int ) ( ( float ) target_total * probgyp );
1880  target_anhydrite = ( int ) ( ( float ) target_total * probgyp * probanh );
1881  target_hemi = ( int ) ( ( float ) target_total * probgyp * probhem );
1882  if ( gsphere(numsize, sphnum, sphrad, sphase) == 1 ) { //unsuccessful generation of microstructure due to excessive amount of particles or too dense
1883  delete [] sphrad;
1884  delete [] sphase;
1885  delete [] sphnum;
1886  return ( 1 );
1887  }
1888  }
1889 
1890  delete [] sphrad;
1891  delete [] sphase;
1892  delete [] sphnum;
1893 
1894 
1895  return ( 0 );
1896 }
1897 
1898 
1899 /* Routine to draw a particle during flocculation routine */
1900 /* See routine chksph for definition of parameters */
1901 void CemhydMatStatus :: drawfloc(int xin, int yin, int zin, int radd, int phasein, int phase2)
1902 {
1903  /* Calls: no other routines */
1904  /* Called by: makefloc */
1905  int xp, yp, zp, i, j, k;
1906  float dist, xdist, ydist, zdist, ftmp;
1907 
1908  /* Check all pixels within the digitized sphere volume */
1909  for ( i = xin - radd; ( i <= xin + radd ); i++ ) {
1910  xp = i;
1911  /* use periodic boundary conditions for sphere placement */
1912  if ( xp < 1 ) {
1913  xp += SYSSIZE;
1914  } else if ( xp > SYSSIZE ) {
1915  xp -= SYSSIZE;
1916  }
1917 
1918  ftmp = ( float ) ( i - xin );
1919  xdist = ftmp * ftmp;
1920  for ( j = yin - radd; ( j <= yin + radd ); j++ ) {
1921  yp = j;
1922  /* use periodic boundary conditions for sphere placement */
1923  if ( yp < 1 ) {
1924  yp += SYSSIZE;
1925  } else if ( yp > SYSSIZE ) {
1926  yp -= SYSSIZE;
1927  }
1928 
1929  ftmp = ( float ) ( j - yin );
1930  ydist = ftmp * ftmp;
1931  for ( k = zin - radd; ( k <= zin + radd ); k++ ) {
1932  zp = k;
1933  /* use periodic boundary conditions for sphere placement */
1934  if ( zp < 1 ) {
1935  zp += SYSSIZE;
1936  } else if ( zp > SYSSIZE ) {
1937  zp -= SYSSIZE;
1938  }
1939 
1940  ftmp = ( float ) ( k - zin );
1941  zdist = ftmp * ftmp;
1942 
1943  /* Compute distance from center of sphere to this pixel */
1944  dist = sqrt(xdist + ydist + zdist);
1945  if ( ( dist - 0.5 ) <= ( float ) radd ) {
1946  /* Update both cement and cemreal images */
1947  cement [ xp ] [ yp ] [ zp ] = phasein;
1948  cemreal [ xp ] [ yp ] [ zp ] = phase2;
1949  }
1950  }
1951  }
1952  }
1953 }
1954 
1955 
1956 /* Routine to check particle placement during flocculation */
1957 /* for particle of size radd centered at (xin,yin,zin) */
1958 /* Returns flag indicating if placement is possible */
1959 int CemhydMatStatus :: chkfloc(int xin, int yin, int zin, int radd)
1960 {
1961  /* Calls: no other routines */
1962  /* Called by: makefloc */
1963  int nofits, xp, yp, zp, i, j, k;
1964  float dist, xdist, ydist, zdist, ftmp;
1965 
1966  nofits = 0; /* Flag indicating if placement is possible */
1967 
1968  /* Check all pixels within the digitized sphere volume */
1969  for ( i = xin - radd; ( ( i <= xin + radd ) && ( nofits == 0 ) ); i++ ) {
1970  xp = i;
1971  /* use periodic boundary conditions for sphere placement */
1972  if ( xp < 1 ) {
1973  xp += SYSSIZE;
1974  } else if ( xp > SYSSIZE ) {
1975  xp -= SYSSIZE;
1976  }
1977 
1978  ftmp = ( float ) ( i - xin );
1979  xdist = ftmp * ftmp;
1980  for ( j = yin - radd; ( ( j <= yin + radd ) && ( nofits == 0 ) ); j++ ) {
1981  yp = j;
1982  /* use periodic boundary conditions for sphere placement */
1983  if ( yp < 1 ) {
1984  yp += SYSSIZE;
1985  } else if ( yp > SYSSIZE ) {
1986  yp -= SYSSIZE;
1987  }
1988 
1989  ftmp = ( float ) ( j - yin );
1990  ydist = ftmp * ftmp;
1991  for ( k = zin - radd; ( ( k <= zin + radd ) && ( nofits == 0 ) ); k++ ) {
1992  zp = k;
1993  /* use periodic boundary conditions for sphere placement */
1994  if ( zp < 1 ) {
1995  zp += SYSSIZE;
1996  } else if ( zp > SYSSIZE ) {
1997  zp -= SYSSIZE;
1998  }
1999 
2000  ftmp = ( float ) ( k - zin );
2001  zdist = ftmp * ftmp;
2002 
2003  /* Compute distance from center of sphere to this pixel */
2004  dist = sqrt(xdist + ydist + zdist);
2005  if ( ( dist - 0.5 ) <= ( float ) radd ) {
2006  if ( ( cement [ xp ] [ yp ] [ zp ] != POROSITY ) ) {
2007  /* Record ID of particle hit */
2008  nofits = cement [ xp ] [ yp ] [ zp ];
2009  }
2010  }
2011 
2012  /* Check for overlap with aggregate */
2013  if ( ( fabs( xp - ( ( float ) ( SYSSIZE + 1 ) / 2.0 ) ) ) < ( ( float ) aggsize / 2.0 ) ) {
2014  nofits = AGG;
2015  }
2016  }
2017  }
2018  }
2019 
2020  /* return flag indicating if sphere will fit */
2021  return ( nofits );
2022 }
2023 
2024 
2025 /* routine to perform flocculation of particles */
2027 /* Calls: drawfloc, chkfloc, ran1 */
2028 /* Called by: main program */
2029 {
2030  int partdo, numfloc;
2031  int nstart;
2032  int nleft = 0, ckall;
2033  int xm, ym, zm, moveran;
2034  int xp, yp, zp, rp, clushit, valkeep;
2035  int iclus;
2036  int *cluspart;
2037  struct cluster *parttmp, *partpoint, *partkeep = NULL;
2038 
2039  cluspart = new int [ NPARTC ];
2040 
2041 
2042  nstart = npart; /* Counter of number of flocs remaining */
2043  for ( iclus = 1; iclus <= npart; iclus++ ) {
2044  cluspart [ iclus ] = iclus;
2045  }
2046 
2047  do {
2048 #ifdef PRINTF
2049  printf("Enter number of flocs desired at end of routine (>0) \n");
2050 #endif
2051 #ifdef CMLFILE
2052  F->get_next_line_in_section(0, ( long & )numfloc);
2053 #endif
2054 #ifdef TINYXML
2055  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, numfloc);
2056 #endif
2057  //fscanf(in, "%d",&numfloc);
2058 #ifdef PRINTF
2059  printf("%d\n", numfloc);
2060 #endif
2061  } while ( numfloc <= 0 );
2062 
2063  while ( nstart > numfloc ) {
2064  nleft = 0;
2065 
2066  /* Try to move each cluster in turn */
2067  for ( iclus = 1; iclus <= npart; iclus++ ) {
2068  if ( clust [ iclus ] == NULL ) {
2069  nleft += 1;
2070  } else {
2071  xm = ym = zm = 0;
2072  /* Generate a random move in one of 6 principal directions */
2073  moveran = ( int ) ( 6. * ran1(seed) );
2074  switch ( moveran ) {
2075  case 0:
2076  xm = 1;
2077  break;
2078  case 1:
2079  xm = ( -1 );
2080  break;
2081  case 2:
2082  ym = 1;
2083  break;
2084  case 3:
2085  ym = ( -1 );
2086  break;
2087  case 4:
2088  zm = 1;
2089  break;
2090  case 5:
2091  zm = ( -1 );
2092  break;
2093  default:
2094  break;
2095  }
2096 
2097  /* First erase all particles in cluster */
2098  partpoint = clust [ iclus ];
2099  while ( partpoint != NULL ) {
2100  xp = partpoint->x;
2101  yp = partpoint->y;
2102  zp = partpoint->z;
2103  rp = partpoint->r;
2104  drawfloc(xp, yp, zp, rp, 0, 0);
2105  partpoint = partpoint->nextpart;
2106  }
2107 
2108  ckall = 0;
2109  /* Now try to draw cluster at new location */
2110  partpoint = clust [ iclus ];
2111  while ( ( partpoint != NULL ) && ( ckall == 0 ) ) {
2112  xp = partpoint->x + xm;
2113  yp = partpoint->y + ym;
2114  zp = partpoint->z + zm;
2115  rp = partpoint->r;
2116  ckall = chkfloc(xp, yp, zp, rp);
2117  partpoint = partpoint->nextpart;
2118  }
2119 
2120  if ( ckall == 0 ) {
2121  /* Place cluster particles at new location */
2122  partpoint = clust [ iclus ];
2123  while ( partpoint != NULL ) {
2124  xp = partpoint->x + xm;
2125  yp = partpoint->y + ym;
2126  zp = partpoint->z + zm;
2127  rp = partpoint->r;
2128  valkeep = partpoint->partphase;
2129  partdo = partpoint->partid;
2130  drawfloc(xp, yp, zp, rp, partdo + CEM - 1, valkeep);
2131  /* Update particle location */
2132  partpoint->x = xp;
2133  partpoint->y = yp;
2134  partpoint->z = zp;
2135  partpoint = partpoint->nextpart;
2136  }
2137  } else {
2138  /* A cluster or aggregate was hit */
2139  /* Draw particles at old location */
2140  partpoint = clust [ iclus ];
2141  /* partkeep stores pointer to last particle in list */
2142  while ( partpoint != NULL ) {
2143  xp = partpoint->x;
2144  yp = partpoint->y;
2145  zp = partpoint->z;
2146  rp = partpoint->r;
2147  valkeep = partpoint->partphase;
2148  partdo = partpoint->partid;
2149  drawfloc(xp, yp, zp, rp, partdo + CEM - 1, valkeep);
2150  partkeep = partpoint;
2151  partpoint = partpoint->nextpart;
2152  }
2153 
2154  /* Determine the cluster hit */
2155  if ( ckall != AGG ) {
2156  clushit = cluspart [ ckall - CEM + 1 ];
2157  /* Move all of the particles from cluster clushit to cluster iclus */
2158  parttmp = clust [ clushit ];
2159  /* Attach new cluster to old one */
2160  partkeep->nextpart = parttmp;
2161  while ( parttmp != NULL ) {
2162  cluspart [ parttmp->partid ] = iclus;
2163  /* Relabel all particles added to this cluster */
2164  parttmp->clustid = iclus;
2165  parttmp = parttmp->nextpart;
2166  }
2167 
2168  /* Disengage the cluster that was hit */
2169  clust [ clushit ] = NULL;
2170  nstart -= 1;
2171  }
2172  }
2173  }
2174  }
2175 
2176  printf("Number left was %d but number of clusters is %d \n", nleft, nstart);
2177  }
2178 
2179  /* end of while loop */
2180  clusleft = nleft;
2181 
2182  delete [] cluspart;
2183 }
2184 
2185 
2186 /* routine to assess global phase fractions present in 3-D system */
2188 /* Calls: no other routines */
2189 /* Called by: main program */
2190 {
2191  long int npor, nc2s, ngyp, ncem, nagg, npozz, ninert, nflyash, nanh, nhem, ncaco3, nslag;
2192  int i, j, k, valph;
2193 
2194  /* counters for the various phase fractions */
2195  npor = 0;
2196  ngyp = 0;
2197  ncem = 0;
2198  nagg = 0;
2199  ninert = 0;
2200  nslag = 0;
2201  nc2s = 0;
2202  npozz = 0;
2203  nflyash = 0;
2204  nanh = 0;
2205  nhem = 0;
2206  ncaco3 = 0;
2207 
2208  /* Check all pixels in 3-D microstructure */
2209  for ( i = 1; i <= SYSSIZE; i++ ) {
2210  for ( j = 1; j <= SYSSIZE; j++ ) {
2211  for ( k = 1; k <= SYSSIZE; k++ ) {
2212  valph = cemreal [ i ] [ j ] [ k ];
2213  if ( valph == POROSITY ) {
2214  npor += 1;
2215  } else if ( valph == CEMID ) {
2216  ncem += 1;
2217  } else if ( valph == C2SID ) {
2218  nc2s += 1;
2219  } else if ( valph == GYPID ) {
2220  ngyp += 1;
2221  } else if ( valph == ANHYDRITE ) {
2222  nanh += 1;
2223  } else if ( valph == HEMIHYDRATE ) {
2224  nhem += 1;
2225  } else if ( valph == AGG ) {
2226  nagg += 1;
2227  } else if ( valph == POZZID ) {
2228  npozz += 1;
2229  } else if ( valph == SLAGID ) {
2230  nslag += 1;
2231  } else if ( valph == INERTID ) {
2232  ninert += 1;
2233  } else if ( valph == FLYASH ) {
2234  nflyash += 1;
2235  } else if ( valph == CACO3 ) {
2236  ncaco3 += 1;
2237  }
2238  }
2239  }
2240  }
2241 
2242  /* Output results */
2243  printf("\n Phase counts are: \n");
2244  printf("Porosity= %ld \n", npor);
2245  printf("Cement= %ld \n", ncem);
2246  printf("C2S= %ld \n", nc2s);
2247  printf("Gypsum= %ld \n", ngyp);
2248  printf("Anhydrite= %ld \n", nanh);
2249  printf("Hemihydrate= %ld \n", nhem);
2250  printf("Pozzolan= %ld \n", npozz);
2251  printf("Inert= %ld \n", ninert);
2252  printf("Slag= %ld \n", nslag);
2253  printf("CaCO3= %ld \n", ncaco3);
2254  printf("Fly Ash= %ld \n", nflyash);
2255  printf("Aggregate= %ld \n", nagg);
2256 }
2257 
2258 
2259 /* Routine to measure phase fractions as a function of distance from */
2260 /* aggregate surface*/
2262 /* Calls: no other routines */
2263 /* Called by: main program */
2264 {
2265  int phase [ 40 ], ptot;
2266  int icnt, ixlo, ixhi, iy, iz, phid, idist;
2267  FILE *aggfile;
2268 
2269  /* By default, results are sent to output file called agglist.out */
2270  aggfile = fopen("agglist.out", "w");
2271  printf("Distance Porosity Cement C2S Gypsum Anhydrite Hemihydrate Pozzolan Inert Slag CaCO3 Fly Ash\n");
2272  fprintf(aggfile, "Distance Porosity Cement C2S Gypsum Anhydrite Hemihydrate Pozzolan Inert Slag CaCO3 Fly Ash\n");
2273 
2274  /* Increase distance from aggregate in increments of one */
2275  for ( idist = 1; idist <= ( SYSSIZE - aggsize ) / 2; idist++ ) {
2276  /* Pixel left of aggregate surface */
2277  ixlo = ( ( SYSSIZE - aggsize + 2 ) / 2 ) - idist;
2278  /* Pixel right of aggregate surface */
2279  ixhi = ( ( SYSSIZE + aggsize ) / 2 ) + idist;
2280 
2281  /* Initialize phase counts for this distance */
2282  for ( icnt = 0; icnt < 39; icnt++ ) {
2283  phase [ icnt ] = 0;
2284  }
2285 
2286  ptot = 0;
2287 
2288  /* Check all pixels which are this distance from aggregate surface */
2289  for ( iy = 1; iy <= SYSSIZE; iy++ ) {
2290  for ( iz = 1; iz <= SYSSIZE; iz++ ) {
2291  phid = cemreal [ ixlo ] [ iy ] [ iz ];
2292  ptot += 1;
2293  if ( phid <= FLYASH ) {
2294  phase [ phid ] += 1;
2295  }
2296 
2297  phid = cemreal [ ixhi ] [ iy ] [ iz ];
2298  ptot += 1;
2299  if ( phid <= FLYASH ) {
2300  phase [ phid ] += 1;
2301  }
2302  }
2303  }
2304 
2305  /* Output results for this distance from surface */
2306  printf("%d %d %d %d %d %d %d %d %d %d %d %d\n", idist, phase [ 0 ], phase [ CEMID ], phase [ C2SID ], phase [ GYPID ], phase [ ANHYDRITE ], phase [ HEMIHYDRATE ], phase [ POZZID ], phase [ INERTID ], phase [ SLAGID ], phase [ CACO3 ], phase [ FLYASH ]);
2307  fprintf(aggfile, "%d %d %d %d %d %d %d %d %d %d %d %d\n", idist, phase [ 0 ], phase [ CEMID ], phase [ C2SID ], phase [ GYPID ], phase [ ANHYDRITE ], phase [ HEMIHYDRATE ], phase [ POZZID ], phase [ INERTID ], phase [ SLAGID ], phase [ CACO3 ], phase [ FLYASH ]);
2308  }
2309 
2310  fclose(aggfile);
2311 }
2312 
2313 
2314 /* routine to assess the connectivity (percolation) of a single phase */
2315 /* Two matrices are used here: one for the current burnt locations */
2316 /* the other to store the newly found burnt locations */
2318 /* Calls: no other routines */
2319 /* Called by: main program */
2320 {
2321  long int inew, ntop, nthrough, ncur, nnew, ntot;
2322  int i, j, k, nmatx [ 29000 ], nmaty [ 29000 ], nmatz [ 29000 ];
2323  int xcn, ycn, zcn, npix, x1, y1, z1, igood, nnewx [ 29000 ], nnewy [ 29000 ], nnewz [ 29000 ];
2324  int jnew, icur;
2325 
2326  do {
2327  printf("Enter phase to analyze 0) pores 1) Cement \n");
2328 #ifdef CMLFILE
2329  F->get_next_line_in_section(0, ( long & )npix);
2330 #endif
2331 #ifdef TINYXML
2332  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, npix);
2333 #endif
2334  //fscanf(in, "%d",&npix);
2335  printf("%d \n", npix);
2336  } while ( ( npix != 0 ) && ( npix != 1 ) );
2337 
2338  /* counters for number of pixels of phase accessible from top surface */
2339  /* and number which are part of a percolated pathway */
2340  ntop = 0;
2341  nthrough = 0;
2342 
2343  /* percolation is assessed from top to bottom only */
2344  /* and burning algorithm is periodic in x and y directions */
2345  k = 1;
2346  for ( i = 1; i <= SYSSIZE; i++ ) {
2347  for ( j = 1; j <= SYSSIZE; j++ ) {
2348  ncur = 0;
2349  ntot = 0;
2350  igood = 0; /* Indicates if bottom has been reached */
2351  if ( ( ( cement [ i ] [ j ] [ k ] == npix ) && ( ( cement [ i ] [ j ] [ SYSSIZE ] == npix ) ||
2352  ( cement [ i ] [ j ] [ SYSSIZE ] == ( npix + BURNTG ) ) ) ) ||
2353  ( ( cement [ i ] [ j ] [ SYSSIZE ] >= CEM ) &&
2354  ( cement [ i ] [ j ] [ k ] >= CEM ) && ( cement [ i ] [ j ] [ k ] < BURNTG ) && ( npix == 1 ) ) ) {
2355  /* Start a burn front */
2356  cement [ i ] [ j ] [ k ] += BURNTG;
2357  ntot += 1;
2358  ncur += 1;
2359  /* burn front is stored in matrices nmat* */
2360  /* and nnew* */
2361  nmatx [ ncur ] = i;
2362  nmaty [ ncur ] = j;
2363  nmatz [ ncur ] = 1;
2364  /* Burn as long as new (fuel) pixels are found */
2365  do {
2366  nnew = 0;
2367  for ( inew = 1; inew <= ncur; inew++ ) {
2368  xcn = nmatx [ inew ];
2369  ycn = nmaty [ inew ];
2370  zcn = nmatz [ inew ];
2371 
2372  /* Check all six neighbors */
2373  for ( jnew = 1; jnew <= 6; jnew++ ) {
2374  x1 = xcn;
2375  y1 = ycn;
2376  z1 = zcn;
2377  if ( jnew == 1 ) {
2378  x1 -= 1;
2379  if ( x1 < 1 ) {
2380  x1 += SYSSIZE;
2381  }
2382  } else if ( jnew == 2 ) {
2383  x1 += 1;
2384  if ( x1 > SYSSIZE ) {
2385  x1 -= SYSSIZE;
2386  }
2387  } else if ( jnew == 3 ) {
2388  y1 -= 1;
2389  if ( y1 < 1 ) {
2390  y1 += SYSSIZE;
2391  }
2392  } else if ( jnew == 4 ) {
2393  y1 += 1;
2394  if ( y1 > SYSSIZE ) {
2395  y1 -= SYSSIZE;
2396  }
2397  } else if ( jnew == 5 ) {
2398  z1 -= 1;
2399  } else if ( jnew == 6 ) {
2400  z1 += 1;
2401  }
2402 
2403  /* Nonperiodic in z direction so be sure to remain in the 3-D box */
2404  if ( ( z1 >= 1 ) && ( z1 <= SYSSIZE ) ) {
2405  if ( ( cement [ x1 ] [ y1 ] [ z1 ] == npix ) || ( ( cement [ x1 ] [ y1 ] [ z1 ] >= CEM ) &&
2406  ( cement [ x1 ] [ y1 ] [ z1 ] < BURNTG ) && ( npix == 1 ) ) ) {
2407  ntot += 1;
2408  cement [ x1 ] [ y1 ] [ z1 ] += BURNTG;
2409  nnew += 1;
2410  if ( nnew >= 29000 ) {
2411  printf("error in size of nnew \n");
2412  }
2413 
2414  nnewx [ nnew ] = x1;
2415  nnewy [ nnew ] = y1;
2416  nnewz [ nnew ] = z1;
2417  /* See if bottom of system has been reached */
2418  if ( z1 == SYSSIZE ) {
2419  igood = 1;
2420  }
2421  }
2422  }
2423  }
2424  }
2425 
2426  if ( nnew > 0 ) {
2427  ncur = nnew;
2428  /* update the burn front matrices */
2429  for ( icur = 1; icur <= ncur; icur++ ) {
2430  nmatx [ icur ] = nnewx [ icur ];
2431  nmaty [ icur ] = nnewy [ icur ];
2432  nmatz [ icur ] = nnewz [ icur ];
2433  }
2434  }
2435  } while ( nnew > 0 );
2436 
2437  ntop += ntot;
2438  if ( igood == 1 ) {
2439  nthrough += ntot;
2440  }
2441  }
2442  }
2443  }
2444 
2445  printf("Phase ID= %d \n", npix);
2446  printf("Number accessible from top= %ld \n", ntop);
2447  printf("Number contained in through pathways= %ld \n", nthrough);
2448 
2449  /* return the burnt sites to their original phase values */
2450  for ( i = 1; i <= SYSSIZE; i++ ) {
2451  for ( j = 1; j <= SYSSIZE; j++ ) {
2452  for ( k = 1; k <= SYSSIZE; k++ ) {
2453  if ( cement [ i ] [ j ] [ k ] >= BURNTG ) {
2454  cement [ i ] [ j ] [ k ] -= BURNTG;
2455  }
2456  }
2457  }
2458  }
2459 }
2460 
2461 
2462 /* Routine to output final microstructure to file */
2464 /* Calls: no other routines */
2465 /* Called by: main program */
2466 {
2467  FILE *outfile, *partfile;
2468  char filen [ 80 ], filepart [ 80 ];
2469  int ix, iy, iz, valout;
2470 
2471 #ifdef PRINTF
2472  printf("Enter name of file to save microstructure to \n");
2473 #endif
2474 #ifdef CMLFILE
2475  F->get_next_line_in_section(0, filen);
2476 #endif
2477 #ifdef TINYXML
2478  QueryStringAttributeExt(xmlFile, "Generate_microstructure", countKey++, filen);
2479 #endif
2480  //fscanf(in, "%s",filen);
2481  printf("%s\n", filen);
2482 
2483  outfile = fopen(filen, "w");
2484 
2485 #ifdef PRINTF
2486  printf("Enter name of file to save particle IDs to \n");
2487 #endif
2488 #ifdef CMLFILE
2489  F->get_next_line_in_section(0, filepart);
2490 #endif
2491 #ifdef TINYXML
2492  QueryStringAttributeExt(xmlFile, "Generate_microstructure", countKey++, filepart);
2493 #endif
2494  //fscanf(in, "%s",filepart);
2495 #ifdef PRINTF
2496  printf("%s\n", filepart);
2497 #endif
2498 
2499  partfile = fopen(filepart, "w");
2500 
2501  for ( iz = 1; iz <= SYSSIZE; iz++ ) {
2502  for ( iy = 1; iy <= SYSSIZE; iy++ ) {
2503  for ( ix = 1; ix <= SYSSIZE; ix++ ) {
2504  valout = cemreal [ ix ] [ iy ] [ iz ];
2505  fprintf(outfile, "%1d\n", valout); //img file
2506  valout = cement [ ix ] [ iy ] [ iz ];
2507  if ( valout < 0 ) {
2508  valout = 0;
2509  }
2510 
2511  fprintf(partfile, "%d\n", valout); //id file
2512  }
2513  }
2514  }
2515 
2516  fclose(outfile);
2517  fclose(partfile);
2518 }
2519 
2520 
2522 {
2523  int userc; /* User choice from menu */
2524  int ig, jg, kg;
2525 
2526  n_sulfate = 0;
2527  target_sulfate = 0;
2528  n_total = 0;
2529  target_total = 0;
2530  n_anhydrite = 0;
2531  target_anhydrite = 0;
2532  n_hemi = 0;
2533  target_hemi = 0;
2534 
2535  alloc_long_3D(cement, SYSIZE + 1);
2537 
2538  clust = new cluster * [ NPARTC ];
2539 
2540  /* Initialize volume array */
2541  volpart [ 0 ] = 1;
2542  volpart [ 1 ] = 19;
2543  volpart [ 2 ] = 81;
2544  volpart [ 3 ] = 179;
2545  volpart [ 4 ] = 389;
2546  volpart [ 5 ] = 739;
2547  volpart [ 6 ] = 1189;
2548  volpart [ 7 ] = 1791;
2549  volpart [ 8 ] = 2553;
2550  volpart [ 9 ] = 3695;
2551  volpart [ 10 ] = 4945;
2552  volpart [ 11 ] = 6403;
2553  volpart [ 12 ] = 8217;
2554  volpart [ 13 ] = 10395;
2555  volpart [ 14 ] = 12893;
2556  volpart [ 15 ] = 15515;
2557  volpart [ 16 ] = 18853;
2558  volpart [ 17 ] = 22575;
2559  volpart [ 18 ] = 26745;
2560  volpart [ 19 ] = 31103;
2561  volpart [ 20 ] = 36137;
2562  volpart [ 21 ] = 41851;
2563  volpart [ 22 ] = 47833;
2564  volpart [ 23 ] = 54435;
2565  volpart [ 24 ] = 61565;
2566  volpart [ 25 ] = 69599;
2567  volpart [ 26 ] = 78205;
2568  volpart [ 27 ] = 87271;
2569  volpart [ 28 ] = 97233;
2570  volpart [ 29 ] = 107783;
2571  volpart [ 30 ] = 119009;
2572  volpart [ 31 ] = 131155;
2573  volpart [ 32 ] = 143761;
2574  volpart [ 33 ] = 157563;
2575  volpart [ 34 ] = 172317;
2576  volpart [ 35 ] = 187511;
2577  volpart [ 36 ] = 203965;
2578 
2579  //printf("Enter random number seed value (a negative integer) \n");
2580  //fscanf(in, "%d",&iseed);
2581  //printf("%d \n",iseed);
2582  nseed = iseed;
2583  seed = & nseed;
2584 
2585  /* Initialize counters and system parameters */
2586  npart = 0;
2587  aggsize = 0;
2588  clusleft = 0;
2589 
2590  /* clear the 3-D system to all porosity to start */
2591  for ( ig = 1; ig <= SYSSIZE; ig++ ) {
2592  for ( jg = 1; jg <= SYSSIZE; jg++ ) {
2593  for ( kg = 1; kg <= SYSSIZE; kg++ ) {
2594  cement [ ig ] [ jg ] [ kg ] = POROSITY; //particle ID
2595  cemreal [ ig ] [ jg ] [ kg ] = POROSITY; //microstructure
2596  }
2597  }
2598  }
2599 
2600  /* present menu and execute user choice */
2601  do {
2602 #ifdef PRINTF
2603  printf(" \n Input User Choice \n");
2604  printf("1) Exit \n");
2605  printf("2) Add spherical particles (cement,gypsum, pozzolans, etc.) to microstructure \n");
2606  printf("3) Flocculate system by reducing number of particle clusters \n");
2607  printf("4) Measure global phase fractions \n");
2608  printf("5) Add an aggregate to the microstructure \n");
2609  printf("6) Measure single phase connectivity (pores or solids) \n");
2610  printf("7) Measure phase fractions vs. distance from aggregate surface \n");
2611  printf("8) Output current microstructure to file \n");
2612 #endif
2613 #ifdef CMLFILE
2614  F->get_next_line_in_section(0, ( long & )userc);
2615 #endif
2616 #ifdef TINYXML
2617  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, userc);
2618 #endif
2619  //fscanf(in, "%d",&userc);
2620  //printf("%d \n",userc);
2621  fflush(stdout);
2622 
2623  switch ( userc ) {
2624  case 2:
2625  if ( create() == 1 ) { //unsuccessful generation of microstructure due to excessive amount of particles or too dense
2626  delete [] clust;
2629  return ( 1 );
2630  }
2631 
2632  break;
2633  case 3:
2634  makefloc();
2635  break;
2636  case 4:
2637  measure();
2638  break;
2639  case 5:
2640  addagg();
2641  break;
2642  case 6:
2643  connect();
2644  break;
2645  case 7:
2646  if ( aggsize != 0 ) {
2647  measagg();
2648  } else {
2649  printf("No aggregate present. \n");
2650  }
2651 
2652  break;
2653  case 8:
2654  outmic();
2655  break;
2656  default:
2657  break;
2658  }
2659  } while ( userc != 1 );
2660 
2661  //store ID in an array
2662  for ( ig = 0; ig < SYSSIZE; ig++ ) {
2663  for ( jg = 0; jg < SYSSIZE; jg++ ) {
2664  for ( kg = 0; kg < SYSSIZE; kg++ ) {
2665  micpart [ ig ] [ jg ] [ kg ] = cement [ ig + 1 ] [ jg + 1 ] [ kg + 1 ];
2666  }
2667  }
2668  }
2669 
2671  delete [] clust;
2672  return ( 0 );
2673 }
2674 
2675 
2676 /*************************************************************************/
2677 /****************************DISTRIB3D************************************/
2678 /*************************************************************************/
2679 
2680 /* routine to create a template for the sphere of interest of radius size */
2681 /* to be used in curvature evaluation */
2682 /* Called by: runsint */
2683 /* Calls no other routines */
2685 {
2686  int icirc, xval, yval, zval;
2687  float xtmp, ytmp;
2688  float dist;
2689 
2690  /* determine and store the locations of all pixels in the 3-D sphere */
2691  icirc = 0;
2692  for ( xval = ( -size ); xval <= size; xval++ ) {
2693  xtmp = ( float ) ( xval * xval );
2694  for ( yval = ( -size ); yval <= size; yval++ ) {
2695  ytmp = ( float ) ( yval * yval );
2696  for ( zval = ( -size ); zval <= size; zval++ ) {
2697  dist = sqrt( xtmp + ytmp + ( float ) ( zval * zval ) );
2698  if ( dist <= ( ( float ) size + 0.5 ) ) {
2699  icirc += 1;
2700  if ( icirc >= MAXSPH ) {
2701  printf("Too many elements in sphere \n");
2702  printf("Must change value of MAXSPH parameter \n");
2703  printf("Currently set at %ld \n", MAXSPH);
2704  exit(1);
2705  }
2706 
2707  xsph [ icirc ] = xval;
2708  ysph [ icirc ] = yval;
2709  zsph [ icirc ] = zval;
2710  }
2711  }
2712  }
2713  }
2714 
2715  /* return the number of pixels contained in sphere of radius (size+0.5) */
2716  return ( icirc );
2717 }
2718 
2719 /* routine to count phase fractions (porosity and solids) */
2720 /* Called by main routine */
2721 /* Calls no other routines */
2723 {
2724  long int npore, nsolid [ 37 ];
2725  int ix, iy, iz;
2726 
2727  npore = 0;
2728  for ( ix = 1; ix < 37; ix++ ) {
2729  nsolid [ ix ] = 0;
2730  }
2731 
2732  /* check all pixels in the 3-D system */
2733  for ( ix = 1; ix <= SYSSIZE; ix++ ) {
2734  for ( iy = 1; iy <= SYSSIZE; iy++ ) {
2735  for ( iz = 1; iz <= SYSSIZE; iz++ ) {
2736  if ( mask [ ix ] [ iy ] [ iz ] == 0 ) {
2737  npore += 1;
2738  } else {
2739  nsolid [ mask [ ix ] [ iy ] [ iz ] ] += 1;
2740  }
2741  }
2742  }
2743  }
2744 
2745  printf("Pores are: %ld \n", npore);
2746  printf("Solids are: %ld %ld %ld %ld %ld %ld\n", nsolid [ 1 ], nsolid [ 2 ],
2747  nsolid [ 3 ], nsolid [ 4 ], nsolid [ 5 ], nsolid [ 6 ]);
2748 }
2749 
2750 /* routine to return number of surface faces exposed to porosity */
2751 /* for pixel located at (xin,yin,zin) */
2752 /* Called by rhcalc */
2753 /* Calls no other routines */
2754 int CemhydMatStatus :: surfpix(int xin, int yin, int zin)
2755 {
2756  int npix, ix1, iy1, iz1;
2757 
2758  npix = 0;
2759 
2760  /* check each of the six immediate neighbors */
2761  /* using periodic boundary conditions */
2762  ix1 = xin - 1;
2763  if ( ix1 < 1 ) {
2764  ix1 += SYSSIZE;
2765  }
2766 
2767  if ( mask [ ix1 ] [ yin ] [ zin ] == 0 ) {
2768  npix += 1;
2769  }
2770 
2771  ix1 = xin + 1;
2772  if ( ix1 > SYSSIZE ) {
2773  ix1 -= SYSSIZE;
2774  }
2775 
2776  if ( mask [ ix1 ] [ yin ] [ zin ] == 0 ) {
2777  npix += 1;
2778  }
2779 
2780  iy1 = yin - 1;
2781  if ( iy1 < 1 ) {
2782  iy1 += SYSSIZE;
2783  }
2784 
2785  if ( mask [ xin ] [ iy1 ] [ zin ] == 0 ) {
2786  npix += 1;
2787  }
2788 
2789  iy1 = yin + 1;
2790  if ( iy1 > SYSSIZE ) {
2791  iy1 -= SYSSIZE;
2792  }
2793 
2794  if ( mask [ xin ] [ iy1 ] [ zin ] == 0 ) {
2795  npix += 1;
2796  }
2797 
2798  iz1 = zin - 1;
2799  if ( iz1 < 1 ) {
2800  iz1 += SYSSIZE;
2801  }
2802 
2803  if ( mask [ xin ] [ yin ] [ iz1 ] == 0 ) {
2804  npix += 1;
2805  }
2806 
2807  iz1 = zin + 1;
2808  if ( iz1 > SYSSIZE ) {
2809  iz1 -= SYSSIZE;
2810  }
2811 
2812  if ( mask [ xin ] [ yin ] [ iz1 ] == 0 ) {
2813  npix += 1;
2814  }
2815 
2816  return ( npix );
2817 }
2818 
2819 /* routine to return the current hydraulic radius for phase phin */
2820 /* Calls surfpix */
2821 /* Called by runsint */
2823 {
2824  int ix, iy, iz;
2825  long int porc, surfc;
2826  float rhval;
2827 
2828  porc = surfc = 0;
2829 
2830  /* Check all pixels in the 3-D volume */
2831  for ( ix = 1; ix <= SYSSIZE; ix++ ) {
2832  for ( iy = 1; iy <= SYSSIZE; iy++ ) {
2833  for ( iz = 1; iz <= SYSSIZE; iz++ ) {
2834  if ( mask [ ix ] [ iy ] [ iz ] == phin ) {
2835  porc += 1;
2836  surfc += surfpix(ix, iy, iz);
2837  }
2838  }
2839  }
2840  }
2841 
2842  printf("Phase area count is %ld \n", porc);
2843  printf("Phase surface count is %ld \n", surfc);
2844  rhval = ( float ) porc * 6. / ( 4. * ( float ) surfc );
2845  printf("Hydraulic radius is %f \n", rhval);
2846  return ( rhval );
2847 }
2848 
2849 /* routine to return count of pixels in a spherical template which are phase */
2850 /* phin or porosity (phase=0) */
2851 /* Calls no other routines */
2852 /* Called by sysinit */
2853 int CemhydMatStatus :: countem(int xp, int yp, int zp, int phin)
2854 {
2855  int xc, yc, zc;
2856  int cumnum, ic;
2857 
2858  cumnum = 0;
2859  for ( ic = 1; ic <= nsph; ic++ ) {
2860  xc = xp + xsph [ ic ];
2861  yc = yp + ysph [ ic ];
2862  zc = zp + zsph [ ic ];
2863  /* Use periodic boundaries */
2864  if ( xc < 1 ) {
2865  xc += SYSSIZE;
2866  } else if ( xc > SYSSIZE ) {
2867  xc -= SYSSIZE;
2868  }
2869 
2870  if ( yc < 1 ) {
2871  yc += SYSSIZE;
2872  } else if ( yc > SYSSIZE ) {
2873  yc -= SYSSIZE;
2874  }
2875 
2876  if ( zc < 1 ) {
2877  zc += SYSSIZE;
2878  } else if ( zc > SYSSIZE ) {
2879  zc -= SYSSIZE;
2880  }
2881 
2882  if ( ( xc != xp ) || ( yc != yp ) || ( zc != zp ) ) {
2883  if ( ( mask [ xc ] [ yc ] [ zc ] == phin ) || ( mask [ xc ] [ yc ] [ zc ] == 0 ) ) {
2884  cumnum += 1;
2885  }
2886  }
2887  }
2888 
2889  return ( cumnum );
2890 }
2891 
2892 /* routine to initialize system by determining local curvature */
2893 /* of all phase 1 and phase 2 pixels */
2894 /* Calls countem */
2895 /* Called by runsint */
2896 void CemhydMatStatus :: sysinit(int ph1, int ph2)
2897 {
2898  int count, xl, yl, zl;
2899 
2900  count = 0;
2901  /* process all pixels in the 3-D box */
2902  for ( xl = 1; xl <= SYSSIZE; xl++ ) {
2903  for ( yl = 1; yl <= SYSSIZE; yl++ ) {
2904  for ( zl = 1; zl <= SYSSIZE; zl++ ) {
2905  /* determine local curvature */
2906  /* For phase 1 want to determine number of porosity pixels */
2907  /* (phase=0) in immediate neighborhood */
2908  if ( mask [ xl ] [ yl ] [ zl ] == ph1 ) {
2909  count = countem(xl, yl, zl, 0);
2910  }
2911 
2912  /* For phase 2 want to determine number of porosity or phase */
2913  /* 2 pixels in immediate neighborhood */
2914  if ( mask [ xl ] [ yl ] [ zl ] == ph2 ) {
2915  count = countem(xl, yl, zl, ph2);
2916  }
2917 
2918  if ( ( count < 0 ) || ( count >= nsph ) ) {
2919  printf("Error count is %d \n", count);
2920  printf("xl %d yl %d zl %d \n", xl, yl, zl);
2921  }
2922 
2923  /* case where we have a phase 1 surface pixel */
2924  /* with non-zero local curvature */
2925  if ( ( count >= 0 ) && ( mask [ xl ] [ yl ] [ zl ] == ph1 ) ) {
2926  curvature [ xl ] [ yl ] [ zl ] = count;
2927  /* update solid curvature histogram */
2928  nsolid [ count ] += 1;
2929  }
2930 
2931  /* case where we have a phase 2 surface pixel */
2932  if ( ( count >= 0 ) && ( mask [ xl ] [ yl ] [ zl ] == ph2 ) ) {
2933  curvature [ xl ] [ yl ] [ zl ] = count;
2934  /* update air curvature histogram */
2935  nair [ count ] += 1;
2936  }
2937  }
2938  }
2939  }
2940 
2941  /* end of xl loop */
2942 }
2943 
2944 /* routine to scan system and determine nsolid (ph2) and nair (ph1) */
2945 /* histograms based on values in phase and curvature arrays */
2946 /* Calls no other routines */
2947 /* Called by runsint */
2948 void CemhydMatStatus :: sysscan(int ph1, int ph2)
2949 {
2950  int xd, yd, zd, curvval;
2951 
2952  /* Scan all pixels in 3-D system */
2953  for ( xd = 1; xd <= SYSSIZE; xd++ ) {
2954  for ( yd = 1; yd <= SYSSIZE; yd++ ) {
2955  for ( zd = 1; zd <= SYSSIZE; zd++ ) {
2956  curvval = curvature [ xd ] [ yd ] [ zd ];
2957 
2958  if ( mask [ xd ] [ yd ] [ zd ] == ph2 ) {
2959  nair [ curvval ] += 1;
2960  } else if ( mask [ xd ] [ yd ] [ zd ] == ph1 ) {
2961  nsolid [ curvval ] += 1;
2962  }
2963  }
2964  }
2965  }
2966 }
2967 
2968 /* routine to return how many cells of solid curvature histogram to use */
2969 /* to accomodate nsearch pixels moving */
2970 /* want to use highest values first */
2971 /* Calls no other routines */
2972 /* Called by movepix */
2974 {
2975  int valfound, i, stop;
2976  long int nsofar;
2977 
2978  /* search histogram from top down until cumulative count */
2979  /* exceeds nsearch */
2980  valfound = nsph - 1;
2981  nsofar = 0;
2982  stop = 0;
2983  for ( i = ( nsph - 1 ); ( ( i >= 0 ) && ( stop == 0 ) ); i-- ) {
2984  nsofar += nsolid [ i ];
2985  if ( nsofar > nsearch ) {
2986  valfound = i;
2987  stop = 1;
2988  }
2989  }
2990 
2991  return ( valfound );
2992 }
2993 
2994 /* routine to determine how many cells of air curvature histogram to use */
2995 /* to accomodate nsearch moving pixels */
2996 /* want to use lowest values first */
2997 /* Calls no other routines */
2998 /* Called by movepix */
2999 
3001 {
3002  int valfound, i, stop;
3003  long int nsofar;
3004 
3005  /* search histogram from bottom up until cumulative count */
3006  /* exceeds nsearch */
3007  valfound = 0;
3008  nsofar = 0;
3009  stop = 0;
3010  for ( i = 0; ( ( i < nsph ) && ( stop == 0 ) ); i++ ) {
3011  nsofar += nair [ i ];
3012  if ( nsofar > nsearch ) {
3013  valfound = i;
3014  stop = 1;
3015  }
3016  }
3017 
3018  return ( valfound );
3019 }
3020 
3021 /* routine to move requested number of pixels (ntomove) from highest */
3022 /* curvature phase 1 (ph1) sites to lowest curvature phase 2 (ph2) sites */
3023 /* Calls procsol and procair */
3024 /* Called by runsint */
3025 int CemhydMatStatus :: movepix(int ntomove, int ph1, int ph2)
3026 {
3027  int count1, count2, ntot, countc, i, xp, yp, zp;
3028  int cmin, cmax, cfg;
3029  int alldone;
3030  long int nsolc, nairc, nsum, nsolm, nairm, nst1, nst2, next1, next2;
3031  float pck, plsol, plair;
3032 
3033  alldone = 0;
3034  /* determine critical values for removal and placement */
3035  count1 = procsol(ntomove);
3036  nsum = 0;
3037  cfg = 0;
3038  cmax = count1;
3039  for ( i = nsph; i > count1; i-- ) {
3040  if ( ( nsolid [ i ] > 0 ) && ( cfg == 0 ) ) {
3041  cfg = 1;
3042  cmax = i;
3043  }
3044 
3045  nsum += nsolid [ i ];
3046  }
3047 
3048  /* Determine movement probability for last cell */
3049  plsol = ( float ) ( ntomove - nsum ) / ( float ) nsolid [ count1 ];
3050  next1 = ntomove - nsum;
3051  nst1 = nsolid [ count1 ];
3052 
3053  count2 = procair(ntomove);
3054  nsum = 0;
3055  cmin = count2;
3056  cfg = 0;
3057  for ( i = 0; i < count2; i++ ) {
3058  if ( ( nair [ i ] > 0 ) && ( cfg == 0 ) ) {
3059  cfg = 1;
3060  cmin = i;
3061  }
3062 
3063  nsum += nair [ i ];
3064  }
3065 
3066  /* Determine movement probability for last cell */
3067  plair = ( float ) ( ntomove - nsum ) / ( float ) nair [ count2 ];
3068  next2 = ntomove - nsum;
3069  nst2 = nair [ count2 ];
3070 
3071  /* Check to see if equilibrium has been reached --- */
3072  /* no further increase in hydraulic radius is possible */
3073  if ( cmin >= cmax ) {
3074  alldone = 1;
3075  printf("Stopping - at equilibrium \n");
3076  printf("cmin- %d cmax- %d \n", cmin, cmax);
3077  return ( alldone );
3078  }
3079 
3080  /* initialize counters for performing sintering */
3081  ntot = 0;
3082  nsolc = 0;
3083  nairc = 0;
3084  nsolm = 0;
3085  nairm = 0;
3086 
3087  /* Now process each pixel in turn */
3088  for ( xp = 1; xp <= SYSSIZE; xp++ ) {
3089  for ( yp = 1; yp <= SYSSIZE; yp++ ) {
3090  for ( zp = 1; zp <= SYSSIZE; zp++ ) {
3091  countc = curvature [ xp ] [ yp ] [ zp ];
3092  /* handle phase 1 case first */
3093  if ( mask [ xp ] [ yp ] [ zp ] == ph1 ) {
3094  if ( countc > count1 ) {
3095  /* convert from phase 1 to phase 2 */
3096  mask [ xp ] [ yp ] [ zp ] = ph2;
3097 
3098  /* update appropriate histogram cells */
3099  nsolid [ countc ] -= 1;
3100  nair [ countc ] += 1;
3101  /* store the location of the modified pixel */
3102  ntot += 1;
3103  }
3104 
3105  if ( countc == count1 ) {
3106  nsolm += 1;
3107  /* generate probability for pixel being removed */
3108  pck = ran1(seed);
3109  if ( ( pck < 0 ) || ( pck > 1.0 ) ) {
3110  pck = 1.0;
3111  }
3112 
3113  if ( ( ( pck < plsol ) && ( nsolc < next1 ) ) || ( ( nst1 - nsolm ) < ( next1 - nsolc ) ) ) {
3114  nsolc += 1;
3115  /* convert phase 1 pixel to phase 2 */
3116  mask [ xp ] [ yp ] [ zp ] = ph2;
3117 
3118  /* update appropriate histogram cells */
3119  nsolid [ count1 ] -= 1;
3120  nair [ count1 ] += 1;
3121  /* store the location of the modified pixel */
3122  ntot += 1;
3123  }
3124  }
3125  }
3126  /* handle phase 2 case here */
3127  else if ( mask [ xp ] [ yp ] [ zp ] == ph2 ) {
3128  if ( countc < count2 ) {
3129  /* convert phase 2 pixel to phase 1 */
3130  mask [ xp ] [ yp ] [ zp ] = ph1;
3131 
3132  nsolid [ countc ] += 1;
3133  nair [ countc ] -= 1;
3134  ntot += 1;
3135  }
3136 
3137  if ( countc == count2 ) {
3138  nairm += 1;
3139  pck = ran1(seed);
3140  if ( ( pck < 0 ) || ( pck > 1.0 ) ) {
3141  pck = 1.0;
3142  }
3143 
3144  if ( ( ( pck < plair ) && ( nairc < next2 ) ) || ( ( nst2 - nairm ) < ( next2 - nairc ) ) ) {
3145  nairc += 1;
3146  /* convert phase 2 to phase 1 */
3147  mask [ xp ] [ yp ] [ zp ] = ph1;
3148 
3149  nsolid [ count2 ] += 1;
3150  nair [ count2 ] -= 1;
3151  ntot += 1;
3152  }
3153  }
3154  }
3155  }
3156 
3157  /* end of zp loop */
3158  }
3159 
3160  /* end of yp loop */
3161  }
3162 
3163  /* end of xloop */
3164  printf("ntot is %d \n", ntot);
3165  return ( alldone );
3166 }
3167 
3168 /* routine to execute user input number of cycles of sintering algorithm */
3169 /* Calls maketemp, rhcalc, sysinit, sysscan, and movepix */
3170 /* Called by main routine */
3171 void CemhydMatStatus :: sinter3d(int ph1id, int ph2id, float rhtarget)
3172 {
3173  int natonce, i, rade, j, rflag;
3174  int keepgo;
3175  long int curvsum1, curvsum2, pixsum1, pixsum2;
3176  float rhnow, avecurv1, avecurv2;
3177 
3178  /* initialize the solid and air count histograms */
3179  for ( i = 0; i <= 499; i++ ) {
3180  nsolid [ i ] = 0;
3181  nair [ i ] = 0;
3182  }
3183 
3184  /* Obtain needed user input */
3185  natonce = 200;
3186  rade = 3;
3187  rflag = 0; /* always initialize system */
3188 
3189  nsph = maketemp(rade);
3190  printf("nsph is %d \n", nsph);
3191  if ( rflag == 0 ) {
3192  sysinit(ph1id, ph2id);
3193  } else {
3194  sysscan(ph1id, ph2id);
3195  }
3196 
3197  i = 0;
3198  rhnow = rhcalc(ph1id);
3199  while ( ( rhnow < rhtarget ) && ( i < MAXCYC_SEAL ) ) {
3200  printf("Now: %f Target: %f \n", rhnow, rhtarget);
3201  i += 1;
3202 #ifdef PRINTF
3203  printf("Cycle: %d \n", i);
3204 #endif
3205  keepgo = movepix(natonce, ph1id, ph2id);
3206  /* If equilibrium is reached, then return to calling routine */
3207  if ( keepgo == 1 ) {
3208  return;
3209  }
3210 
3211  curvsum1 = 0;
3212  curvsum2 = 0;
3213  pixsum1 = 0;
3214  pixsum2 = 0;
3215  /* Determine average curvatures for phases 1 and 2 */
3216  for ( j = 0; j <= nsph; j++ ) {
3217  pixsum1 += nsolid [ j ];
3218  curvsum1 += ( j * nsolid [ j ] );
3219  pixsum2 += nair [ j ];
3220  curvsum2 += ( j * nair [ j ] );
3221  }
3222 
3223  avecurv1 = ( float ) curvsum1 / ( float ) pixsum1;
3224  avecurv2 = ( float ) curvsum2 / ( float ) pixsum2;
3225  printf("Ave. solid curvature: %f \n", avecurv1);
3226  printf("Ave. air curvature: %f \n", avecurv2);
3227  rhnow = rhcalc(ph1id);
3228  }
3229 }
3230 
3232 {
3233  int valin, ix, iy, iz;
3234  int ix1, iy1, iz1, k;
3235  long int voltot, surftot;
3236 
3237  for ( ix = 0; ix <= 42; ix++ ) {
3238  volume [ ix ] = surface [ ix ] = 0;
3239  }
3240 
3241  /* Read in image and accumulate volume totals */
3242  for ( iz = 1; iz <= SYSIZE; iz++ ) {
3243  for ( iy = 1; iy <= SYSIZE; iy++ ) {
3244  for ( ix = 1; ix <= SYSIZE; ix++ ) {
3245  valin = mask [ ix ] [ iy ] [ iz ];
3246  volume [ valin ] += 1;
3247  }
3248  }
3249  }
3250 
3251 
3252  for ( iz = 1; iz <= SYSIZE; iz++ ) {
3253  for ( iy = 1; iy <= SYSIZE; iy++ ) {
3254  for ( ix = 1; ix <= SYSIZE; ix++ ) {
3255  if ( mask [ ix ] [ iy ] [ iz ] != 0 ) {
3256  valin = mask [ ix ] [ iy ] [ iz ];
3257  /* Check six neighboring pixels for porosity */
3258  for ( k = 1; k <= 6; k++ ) {
3259  switch ( k ) {
3260  case 1:
3261  ix1 = ix - 1;
3262  if ( ix1 < 1 ) {
3263  ix1 += SYSIZE;
3264  }
3265 
3266  iy1 = iy;
3267  iz1 = iz;
3268  break;
3269  case 2:
3270  ix1 = ix + 1;
3271  if ( ix1 > SYSIZE ) {
3272  ix1 -= SYSIZE;
3273  }
3274 
3275  iy1 = iy;
3276  iz1 = iz;
3277  break;
3278  case 3:
3279  iy1 = iy - 1;
3280  if ( iy1 < 1 ) {
3281  iy1 += SYSIZE;
3282  }
3283 
3284  ix1 = ix;
3285  iz1 = iz;
3286  break;
3287  case 4:
3288  iy1 = iy + 1;
3289  if ( iy1 > SYSIZE ) {
3290  iy1 -= SYSIZE;
3291  }
3292 
3293  ix1 = ix;
3294  iz1 = iz;
3295  break;
3296  case 5:
3297  iz1 = iz - 1;
3298  if ( iz1 < 1 ) {
3299  iz1 += SYSIZE;
3300  }
3301 
3302  iy1 = iy;
3303  ix1 = ix;
3304  break;
3305  case 6:
3306  iz1 = iz + 1;
3307  if ( iz1 > SYSIZE ) {
3308  iz1 -= SYSIZE;
3309  }
3310 
3311  iy1 = iy;
3312  ix1 = ix;
3313  break;
3314  default:
3315  break;
3316  }
3317 
3318  if ( ( ix1 < 1 ) || ( iy1 < 1 ) || ( iz1 < 1 ) || ( ix1 > SYSIZE ) || ( iy1 > SYSIZE ) || ( iz1 > SYSIZE ) ) {
3319  printf("%d %d %d \n", ix1, iy1, iz1);
3320  exit(1);
3321  }
3322 
3323  if ( mask [ ix1 ] [ iy1 ] [ iz1 ] == 0 ) {
3324  surface [ valin ] += 1;
3325  }
3326  }
3327  }
3328  }
3329  }
3330  }
3331 
3332 #ifdef PRINTF
3333  printf("Phase Volume Surface Volume Surface \n");
3334  printf(" ID count count fraction fraction \n");
3335 #endif
3336  /* Only include clinker phases in surface area fraction calculation */
3337  surftot = surface [ 1 ] + surface [ 2 ] + surface [ 3 ] + surface [ 4 ];
3338  voltot = volume [ 1 ] + volume [ 2 ] + volume [ 3 ] + volume [ 4 ];
3339  k = 0;
3340 #ifdef PRINTF
3341  printf(" %d %8ld %8ld \n", k, volume [ 0 ], surface [ 0 ]);
3342 
3343  for ( k = 1; k <= 4; k++ ) {
3344  printf(" %d %8ld %8ld %.5f %.5f\n", k, volume [ k ], surface [ k ],
3345  ( float ) volume [ k ] / ( float ) voltot, ( float ) surface [ k ] / ( float ) surftot);
3346  }
3347 
3348  printf("Total %8ld %8ld\n\n\n", voltot, surftot);
3349 
3350  for ( k = 5; k <= 11; k++ ) {
3351  printf(" %d %8ld %8ld\n", k, volume [ k ], surface [ k ]);
3352  }
3353 
3354  printf(" 20 %8ld %8ld\n", volume [ 20 ], surface [ 20 ]);
3355 
3356  for ( k = 24; k <= 27; k++ ) {
3357  printf(" %d %8ld %8ld\n", k, volume [ k ], surface [ k ]);
3358  }
3359 
3360  printf(" 28 %8ld %8ld\n", volume [ 28 ], surface [ 28 ]);
3361 #endif
3362  ( void ) surftot;
3363  ( void ) voltot;
3364 }
3365 
3366 void CemhydMatStatus :: rand3d(int phasein, int phaseout, float xpt)
3367 {
3368  float s2, ss, sdiff, xtmp, ytmp;
3369  //static float normm[SYSIZE+1][SYSIZE+1][SYSIZE+1];
3370  //static float res[SYSIZE+1][SYSIZE+1][SYSIZE+1];
3371  double ***normm, ***res;
3372  //static float filter [32][32][32];
3373  double ***filter;
3374  int done, r [ 61 ];
3375  //static float s[61],xr[61],sum[502];
3376  float *s, *xr, *sum;
3377  double val2;
3378  double t1, t2, x1, x2, u1, u2, xrad, resmax, resmin;
3379  float xtot, filval, radius, sect, sumtot, vcrit;
3380  int valin, r1, r2, i1, i2, i3, i, j, k, j1, k1;
3381  int ido, iii, jjj, ix, iy, iz, index;
3382 #ifdef CMLFILE
3383  char tempstr [ 256 ];
3384 #endif
3385 
3386  alloc_double_3D(normm, SYSIZE + 1);
3387  alloc_double_3D(res, SYSIZE + 1);
3388  alloc_double_3D(filter, 32);
3389 
3390  s = new float [ 61 ];
3391  xr = new float [ 61 ];
3392  sum = new float [ 502 ];
3393 
3394  /* Create the Gaussian noise image */
3395  i1 = i2 = i3 = 1;
3396  for ( i = 1; i <= ( ( SYSIZE * SYSIZE * SYSIZE ) / 2 ); i++ ) {
3397  u1 = ran1(seed);
3398  u2 = ran1(seed);
3399  t1 = 2. * M_PI * u2;
3400  t2 = sqrt( -2. * log(u1) );
3401  x1 = cos(t1) * t2;
3402  x2 = sin(t1) * t2;
3403  normm [ i1 ] [ i2 ] [ i3 ] = x1;
3404  i1 += 1;
3405  if ( i1 > SYSIZE ) {
3406  i1 = 1;
3407  i2 += 1;
3408  if ( i2 > SYSIZE ) {
3409  i2 = 1;
3410  i3 += 1;
3411  }
3412  }
3413 
3414  normm [ i1 ] [ i2 ] [ i3 ] = x2;
3415  i1 += 1;
3416  if ( i1 > SYSIZE ) {
3417  i1 = 1;
3418  i2 += 1;
3419  if ( i2 > SYSIZE ) {
3420  i2 = 1;
3421  i3 += 1;
3422  }
3423  }
3424  }
3425 
3426  /* Now perform the convolution */
3427 #ifdef CMLFILE
3428  F->get_next_line_in_section(0, ( long & )ido);
3429 #endif
3430  //fscanf(in,"%d",&ido);
3431 #ifdef TINYXML
3432  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, ido);
3433 #endif
3434 
3435 #ifdef PRINTF
3436  printf("Number of points in correlation file is %d \n", ido);
3437 #endif
3438 
3439  for ( i = 1; i <= ido; i++ ) {
3440 #ifdef CMLFILE
3441  F->get_next_line_in_section(0, tempstr);
3442  sscanf(tempstr, "%d %f", & valin, & val2);
3443 #endif
3444 #ifdef TINYXML
3445  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, valin);
3446  QueryNumAttributeExt(xmlFile, "Generate_microstructure", countKey++, val2);
3447 #endif
3448  //fscanf(in,"%d %f",&valin,&val2);
3449  r [ i ] = valin;
3450  s [ i ] = ( float ) val2;
3451  xr [ i ] = ( float ) r [ i ];
3452  }
3453 
3454  ss = s [ 1 ];
3455  s2 = ss * ss;
3456  /* Load up the convolution matrix */
3457  sdiff = ss - s2;
3458  for ( i = 0; i < 31; i++ ) {
3459  iii = i * i;
3460  for ( j = 0; j < 31; j++ ) {
3461  jjj = j * j;
3462  for ( k = 0; k < 31; k++ ) {
3463  xtmp = ( float ) ( iii + jjj + k * k );
3464  radius = sqrt(xtmp);
3465  r1 = ( int ) radius + 1;
3466  r2 = r1 + 1;
3467  if ( s [ r1 ] < 0.0 ) {
3468  printf("%d and %d %f and %f with xtmp of %f\n", r1, r2, s [ r1 ], s [ r2 ], xtmp);
3469  fflush(stdout);
3470  exit(1);
3471  }
3472 
3473  xrad = radius + 1 - r1;
3474  filval = s [ r1 ] + ( s [ r2 ] - s [ r1 ] ) * xrad;
3475  filter [ i + 1 ] [ j + 1 ] [ k + 1 ] = ( filval - s2 ) / sdiff;
3476  }
3477  }
3478  }
3479 
3480  /* Now filter the image maintaining periodic boundaries */
3481  /*fixed for periodic boundaries (small microstructures) - smilauer 4.12.2006*/
3482  resmax = 0.0;
3483  resmin = 1.0;
3484  for ( i = 1; i <= SYSIZE; i++ ) {
3485  for ( j = 1; j <= SYSIZE; j++ ) {
3486  for ( k = 1; k <= SYSIZE; k++ ) {
3487  res [ i ] [ j ] [ k ] = 0.0;
3488  if ( ( float ) mask [ i ] [ j ] [ k ] == phasein ) {
3489  for ( ix = 1; ix <= 31; ix++ ) {
3490  i1 = i + ix - 1;
3491  while ( i1 < 1 ) { //if(i1<1){i1+=SYSIZE;}
3492  i1 += SYSIZE;
3493  }
3494 
3495  while ( i1 > SYSIZE ) { //else if(i1>SYSIZE){
3496  i1 -= SYSIZE;
3497  }
3498 
3499  for ( iy = 1; iy <= 31; iy++ ) {
3500  j1 = j + iy - 1;
3501  while ( j1 < 1 ) { //if(j1<1){j1+=SYSIZE;}
3502  j1 += SYSIZE;
3503  }
3504 
3505  while ( j1 > SYSIZE ) { //else if(j1>SYSIZE){
3506  j1 -= SYSIZE;
3507  }
3508 
3509  for ( iz = 1; iz <= 31; iz++ ) {
3510  k1 = k + iz - 1;
3511  while ( k1 < 1 ) { //if(k1<1){k1+=SYSIZE;}
3512  k1 += SYSIZE;
3513  }
3514 
3515  while ( k1 > SYSIZE ) { //else if(k1>SYSIZE){
3516  k1 -= SYSIZE;
3517  }
3518 
3519  res [ i ] [ j ] [ k ] += normm [ i1 ] [ j1 ] [ k1 ] * filter [ ix ] [ iy ] [ iz ];
3520  }
3521  }
3522  }
3523 
3524  if ( res [ i ] [ j ] [ k ] > resmax ) {
3525  resmax = res [ i ] [ j ] [ k ];
3526  }
3527 
3528  if ( res [ i ] [ j ] [ k ] < resmin ) {
3529  resmin = res [ i ] [ j ] [ k ];
3530  }
3531  }
3532  }
3533 
3534 #ifdef PRINTF
3535  printf(".");
3536 #endif
3537  }
3538 
3539 #ifdef PRINTF
3540  printf("%d out of %d\n", i, SYSIZE);
3541 #endif
3542  }
3543 
3544 #ifdef PRINTF
3545  printf("\n");
3546 #endif
3547  /* Now threshold the image */
3548  sect = ( resmax - resmin ) / 500.;
3549  for ( i = 1; i <= 500; i++ ) {
3550  sum [ i ] = 0.0;
3551  }
3552 
3553  xtot = 0.0;
3554  for ( i = 1; i <= SYSIZE; i++ ) {
3555  for ( j = 1; j <= SYSIZE; j++ ) {
3556  for ( k = 1; k <= SYSIZE; k++ ) {
3557  if ( ( float ) mask [ i ] [ j ] [ k ] == phasein ) {
3558  xtot += 1.0;
3559  index = 1 + ( int ) ( ( res [ i ] [ j ] [ k ] - resmin ) / sect );
3560  if ( index > 500 ) {
3561  index = 500;
3562  }
3563 
3564  sum [ index ] += 1.0;
3565  }
3566  }
3567  }
3568  }
3569 
3570  /* Determine which bin to choose for correct thresholding */
3571  sumtot = vcrit = 0.0;
3572  done = 0;
3573  for ( i = 1; ( ( i <= 500 ) && ( done == 0 ) ); i++ ) {
3574  sumtot += sum [ i ] / xtot;
3575  if ( sumtot > xpt ) {
3576  ytmp = ( float ) i;
3577  vcrit = resmin + ( resmax - resmin ) * ( ytmp - 0.5 ) / 500.;
3578  done = 1;
3579  }
3580  }
3581 
3582 #ifdef PRINTF
3583  printf("Critical volume fraction is %f\n", vcrit);
3584 #endif
3585 
3586  for ( k = 1; k <= SYSIZE; k++ ) {
3587  for ( j = 1; j <= SYSIZE; j++ ) {
3588  for ( i = 1; i <= SYSIZE; i++ ) {
3589  if ( ( float ) mask [ i ] [ j ] [ k ] == phasein ) {
3590  if ( res [ i ] [ j ] [ k ] > vcrit ) {
3591  mask [ i ] [ j ] [ k ] = phaseout;
3592  }
3593  }
3594  }
3595  }
3596  }
3597 
3598  dealloc_double_3D(normm, SYSIZE + 1);
3599  dealloc_double_3D(res, SYSIZE + 1);
3600  dealloc_double_3D(filter, 32);
3601 
3602  delete [] s;
3603  delete [] xr;
3604  delete [] sum;
3605 }
3606 
3607 
3608 /*disabled sintering*/
3610 {
3611  int i, j, k, alumval, alum2, valin;
3612  int output_img;
3613  double volin, volf [ 5 ], surff [ 5 ], rhtest, rdesire;
3614  char filen [ 80 ];
3615  //char fileout[80],filecem[80]
3616  FILE *infile;
3617  FILE *outfile_img = NULL, *outfile_id = NULL;
3618 
3620 
3621  /* Seed the random number generator */
3622  //printf("Enter random number seed (negative integer) \n");
3623  //fscanf(in, "%d",&nseed);
3624  //printf("%d\n",nseed);
3625  nseed = iseed;
3626 #ifdef PRINTF
3627  printf("%d\n", * seed);
3628 #endif
3629 
3630  /* Read in the parameters to use */
3631  //printf("Enter name of cement microstructure image file\n");
3632  //fscanf(in, "%s",filen);
3633  //printf("%s\n",filen);
3634 
3635 
3636  /* Set up the correlation filenames
3637  * printf("Enter name of sil correlation files\n");
3638  * fscanf(in, "%s",filesil);
3639  * printf("%s\n",filesil);
3640  * printf("Enter name of c3s correlation files\n");
3641  * fscanf(in, "%s",filec3s);
3642  * printf("%s\n",filec3s);
3643  * printf("Enter name of c4af correlation files\n");
3644  * fscanf(in, "%s",filealum);
3645  * printf("%s\n",filealum);
3646  */
3647 
3648  alumval = 4;
3649 
3650  //assume always C4AF
3651  /*
3652  * testfile=fopen(filealum,"r");
3653  * if(testfile==NULL){
3654  * alumflag=0;
3655  * sprintf(filealum,"%s",filecem);
3656  * strcat(filealum,".c3a");
3657  * alumval=3;
3658  * }
3659  * else{
3660  * fclose(testfile);
3661  * }
3662  */
3663  //printf("Enter name of new cement microstructure image file\n");
3664  //fscanf(in, "%s",fileout);
3665  // printf("%s\n",fileout);
3666  for ( i = 1; i <= 4; i++ ) {
3667 #ifdef CMLFILE
3668  F->get_next_line_in_section(0, volin);
3669 #endif
3670 #ifdef TINYXML
3671  if ( i == 1 ) {
3672  QueryNumAttributeExt(xmlFile, "Generate_microstructure", "C3S_unit_frac", volin);
3673  }
3674 
3675  if ( i == 2 ) {
3676  QueryNumAttributeExt(xmlFile, "Generate_microstructure", "C2S_unit_frac", volin);
3677  }
3678 
3679  if ( i == 3 ) {
3680  QueryNumAttributeExt(xmlFile, "Generate_microstructure", "C3A_unit_frac", volin);
3681  }
3682 
3683  if ( i == 4 ) {
3684  QueryNumAttributeExt(xmlFile, "Generate_microstructure", "C4AF_unit_frac", volin);
3685  }
3686 
3687 #endif
3688  //fscanf(in, "%f",&volin);
3689  volf [ i ] = volin;
3690 #ifdef PRINTF
3691  printf("Volume %f\n", volf [ i ]);
3692 #endif
3693  //fscanf(in, "%f",&volin);
3694  surff [ i ] = volin;
3695 #ifdef PRINTF
3696  printf("Surface %f\n", surff [ i ]);
3697 #endif
3698  }
3699 
3700  /* Read in the original microstructure image file */
3701  if ( ( infile = fopen(filen, "r") ) != NULL ) {
3702  for ( k = 1; k <= SYSIZE; k++ ) {
3703  for ( j = 1; j <= SYSIZE; j++ ) {
3704  for ( i = 1; i <= SYSIZE; i++ ) {
3705  if ( fscanf(infile, "%d", & valin) != 1 ) {
3706  OOFEM_ERROR("distrib3d reading error");
3707  }
3708 
3709  mask [ i ] [ j ] [ k ] = valin;
3710  curvature [ i ] [ j ] [ k ] = 0;
3711  }
3712  }
3713  }
3714 
3715  fclose(infile);
3716  } else {
3717  for ( k = 1; k <= SYSIZE; k++ ) {
3718  for ( j = 1; j <= SYSIZE; j++ ) {
3719  for ( i = 1; i <= SYSIZE; i++ ) {
3720  mask [ i ] [ j ] [ k ] = cemreal [ i ] [ j ] [ k ];
3721  }
3722  }
3723  }
3724  }
3725 
3726  stat3d();
3727 
3728  /* First filtering */
3729  volin = volf [ 1 ] + volf [ 2 ];
3730  if ( volin < 1.0 ) {
3731  rand3d(1, alumval, volin);
3732 
3733  /* First sintering */
3734  stat3d();
3735  rdesire = ( surff [ 1 ] + surff [ 2 ] ) * ( float ) ( surface [ 1 ] + surface [ alumval ] );
3736  if ( rdesire != 0.0 ) {
3737  if ( ( int ) rdesire < surface [ 1 ] ) {
3738  rhtest = ( 6. / 4. ) * ( float ) ( volume [ 1 ] ) / rdesire;
3739  //sinter3d(1,alumval,rhtest);
3740  } else {
3741  rdesire = ( surff [ 3 ] + surff [ 4 ] ) * ( float ) ( surface [ 1 ] + surface [ alumval ] );
3742  if ( rdesire != 0.0 ) {
3743  rhtest = ( 6. / 4. ) * ( float ) ( volume [ alumval ] ) / rdesire;
3744  //sinter3d(alumval,1,rhtest);
3745  }
3746  }
3747  }
3748  }
3749 
3750  /* Second filtering */
3751  if ( ( volf [ 1 ] + volf [ 2 ] ) > 0.0 ) {
3752  volin = volf [ 1 ] / ( volf [ 1 ] + volf [ 2 ] );
3753  if ( volin < 1.0 ) {
3754  rand3d(1, 2, volin);
3755 
3756  /* Second sintering */
3757  stat3d();
3758  rdesire = ( surff [ 1 ] / ( surff [ 1 ] + surff [ 2 ] ) ) * ( float ) ( surface [ 1 ] + surface [ 2 ] );
3759  if ( rdesire != 0.0 ) {
3760  if ( ( int ) rdesire < surface [ 1 ] ) {
3761  rhtest = ( 6. / 4. ) * ( float ) ( volume [ 1 ] ) / rdesire;
3762  //sinter3d(1,2,rhtest);
3763  } else {
3764  rdesire = ( surff [ 2 ] / ( surff [ 1 ] + surff [ 2 ] ) ) * ( float ) ( surface [ 1 ] + surface [ 2 ] );
3765  if ( rdesire != 0.0 ) {
3766  rhtest = ( 6. / 4. ) * ( float ) ( volume [ 2 ] ) / rdesire;
3767  //sinter3d(2,1,rhtest);
3768  }
3769  }
3770  }
3771  }
3772  }
3773 
3774  /* Third (final) filtering */
3775  if ( alumval == 4 ) {
3776  volin = volf [ 4 ] / ( volf [ 4 ] + volf [ 3 ] );
3777  alum2 = 3;
3778  } else {
3779  volin = volf [ 3 ] / ( volf [ 4 ] + volf [ 3 ] );
3780  alum2 = 4;
3781  }
3782 
3783  if ( volin < 1.0 ) {
3784  rand3d(alumval, alum2, volin);
3785 
3786  /* Third (final) sintering */
3787  stat3d();
3788  if ( alumval == 4 ) {
3789  rdesire = ( surff [ 4 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) ( surface [ 3 ] + surface [ 4 ] );
3790  if ( rdesire != 0.0 ) {
3791  if ( ( int ) rdesire < surface [ 4 ] ) {
3792  rhtest = ( 6. / 4. ) * ( float ) ( volume [ 4 ] ) / rdesire;
3793  //sinter3d(alumval,alum2,rhtest);
3794  } else {
3795  rdesire = ( surff [ 3 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) ( surface [ 3 ] + surface [ 4 ] );
3796  if ( rdesire != 0.0 ) {
3797  rhtest = ( 6. / 4. ) * ( float ) ( volume [ 3 ] ) / rdesire;
3798  //sinter3d(alum2,alumval,rhtest);
3799  }
3800  }
3801  }
3802  } else {
3803  rdesire = ( surff [ 3 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) ( surface [ 3 ] + surface [ 4 ] );
3804  if ( rdesire != 0.0 ) {
3805  if ( ( int ) rdesire < surface [ 3 ] ) {
3806  rhtest = ( 6. / 4. ) * ( float ) ( volume [ 3 ] ) / rdesire;
3807  //sinter3d(alumval,alum2,rhtest);
3808  } else {
3809  rdesire = ( surff [ 4 ] / ( surff [ 3 ] + surff [ 4 ] ) ) * ( float ) ( surface [ 3 ] + surface [ 4 ] );
3810  if ( rdesire != 0.0 ) {
3811  rhtest = ( 6. / 4. ) * ( float ) ( volume [ 4 ] ) / rdesire;
3812  //sinter3d(alum2,alumval,rhtest);
3813  }
3814  }
3815  }
3816  }
3817  }
3818 
3819  /* Output final microstructure */
3820 
3821  //check whether to save it into img and id files
3822 #ifdef CMLFILE
3823  F->get_value(14, ( long & )output_img);
3824 #endif
3825 #ifdef TINYXML
3826  QueryNumAttributeExt(xmlFile, "Output_initial_microstructure", 0, output_img);
3827 #endif
3828 
3829 
3830  if ( output_img ) {
3831 #ifdef CMLFILE
3832  F->get_value(15, filen);
3833 #endif
3834 #ifdef TINYXML
3835  QueryStringAttributeExt(xmlFile, "Output_initial_microstructure_img_file", 0, filen);
3836 #endif
3837  if ( ( outfile_img = fopen(filen, "w") ) == NULL ) {
3838  printf("Output img file %s can not be created (file %s, line %d)\n", filen, __FILE__, __LINE__);
3839  exit(1);
3840  }
3841 
3842 #ifdef CMLFILE
3843  F->get_value(16, filen);
3844 #endif
3845 #ifdef TINYXML
3846  QueryStringAttributeExt(xmlFile, "Output_initial_microstructure_id_file", 0, filen);
3847 #endif
3848  if ( ( outfile_id = fopen(filen, "w") ) == NULL ) {
3849  printf("Output id file %s can not be created (file %s, line %d)\n", filen, __FILE__, __LINE__);
3850  exit(1);
3851  }
3852  }
3853 
3854  for ( k = 0; k < SYSIZE; k++ ) {
3855  for ( j = 0; j < SYSIZE; j++ ) {
3856  for ( i = 0; i < SYSIZE; i++ ) {
3857  mic [ i ] [ j ] [ k ] = mask [ i + 1 ] [ j + 1 ] [ k + 1 ];
3858  micorig [ i ] [ j ] [ k ] = mic [ i ] [ j ] [ k ];
3859  if ( output_img ) {
3860  fprintf(outfile_img, "%d\n", mic [ i ] [ j ] [ k ]);
3861  fprintf(outfile_id, "%ld\n", micpart [ i ] [ j ] [ k ]);
3862  }
3863  }
3864  }
3865  }
3866 
3867  if ( output_img ) {
3868  fclose(outfile_img);
3869  fclose(outfile_id);
3870  }
3871 
3872  //deallocate curvature
3873  dealloc_int_3D(curvature, SYSIZE + 1);
3874 
3875  //deallocate cemreal[][][]
3876  dealloc_long_3D(cemreal, SYSIZE + 1);
3877  ( void ) rhtest;
3878 }
3879 
3880 
3881 /**********************************************************/
3882 /*************DISREALNEW***********************************/
3883 /**********************************************************/
3884 
3886 {
3887  int i;
3888  double slagin, CHperslag;
3889  FILE *slagfile;
3890 
3891  for ( i = 0; i <= EMPTYP; i++ ) {
3892  creates [ i ] = 0;
3893  soluble [ i ] = 0;
3894  disprob [ i ] = 0.0;
3895  disbase [ i ] = 0.0;
3896  pHeffect [ i ] = 0.0;
3897  }
3898 
3899  /* soluble [x] - flag indicating if phase x is soluble */
3900  /* disprob [x] - probability of dissolution (relative diss. rate) */
3901  gypabsprob = 0.01; /* One sulfate absorbed per 100 CSH units */
3902  /* Source is from Taylor's Cement Chemistry for K and Na absorption */
3903  /* Note that for the first cycle, of the clinker phases only the */
3904  /* aluminates and gypsum are soluble */
3905  soluble [ C4AF ] = 1;
3906  disprob [ C4AF ] = disbase [ C4AF ] = 0.067 / DISBIAS;
3907  creates [ C4AF ] = POROSITY;
3908  pHeffect [ C4AF ] = 1.0;
3909  soluble [ C3S ] = 0;
3910  disprob [ C3S ] = disbase [ C3S ] = 0.7 / DISBIAS;
3911  creates [ C3S ] = DIFFCSH;
3912  pHeffect [ C3S ] = 1.0;
3913  soluble [ C2S ] = 0;
3914  disprob [ C2S ] = disbase [ C2S ] = 0.1 / DISBIAS;
3915  creates [ C2S ] = DIFFCSH;
3916  pHeffect [ C2S ] = 1.0;
3917  soluble [ C3A ] = 1;
3918  /* increased back to 0.4 from 0.25 7/8/99 */
3919  disprob [ C3A ] = disbase [ C3A ] = 0.4 / DISBIAS;
3920  creates [ C3A ] = POROSITY;
3921  pHeffect [ C3A ] = 1.0;
3922  soluble [ GYPSUM ] = 1;
3923  /* Changed from 0.05 to 0.015 9/29/98 */
3924  /* Changed to 0.040 10/15/98 */
3925  /* back to 0.05 from 0.10 7/8/99 */
3926  /* from 0.05 to 0.02 4/4/00 */
3927  /* from 0.02 to 0.025 8/13/01 */
3928  disprob [ GYPSUM ] = disbase [ GYPSUM ] = 0.025; /*geaendert am 04.04.00, urspr. 0.05*/
3929  /* dissolved gypsum distributed at random throughout microstructure */
3930  creates [ GYPSUM ] = POROSITY;
3931  soluble [ GYPSUMS ] = 1;
3932  pHeffect [ GYPSUM ] = 0.0;
3933  /* Changed from 0.05 to 0.015 9/29/98 */
3934  /* Changed to 0.020 10/15/98 */
3935  /* and also changed all sulfate based dissolution rates */
3936  disprob [ GYPSUMS ] = disbase [ GYPSUMS ] = disprob [ GYPSUM ];
3937  creates [ GYPSUMS ] = POROSITY;
3938  pHeffect [ GYPSUMS ] = 0.0;
3939  soluble [ ANHYDRITE ] = 1;
3940  /* set anhydrite dissolution at 4/5ths of that of gypsum */
3941  /* Source: Uchikawa et al., CCR, 1984 */
3942  disprob [ ANHYDRITE ] = disbase [ ANHYDRITE ] = 0.8 * disprob [ GYPSUM ];
3943  /* dissolved anhydrite distributed at random throughout microstructure */
3944  creates [ ANHYDRITE ] = POROSITY;
3945  pHeffect [ ANHYDRITE ] = 0.0;
3946  soluble [ HEMIHYD ] = 1;
3947  /* set hemihydrate dissolution at 3 times that of gypsum */
3948  /* Source: Uchikawa et al., CCR, 1984 */
3949  /* Changed to 1.5 times that of gypsum 6/1/00 */
3950  disprob [ HEMIHYD ] = disbase [ HEMIHYD ] = 1.5 * disprob [ GYPSUM ]; /* ge´┐Żndert am 01.06.00, urspr. 3.0 */
3951  /* dissolved hemihydrate distributed at random throughout microstructure */
3952  creates [ HEMIHYD ] = POROSITY;
3953  pHeffect [ HEMIHYD ] = 0.0;
3954  /* CH soluble to allow for Ostwald ripening of crystals */
3955  soluble [ CH ] = 1;
3956  disprob [ CH ] = disbase [ CH ] = 0.5 / DISBIAS;
3957  creates [ CH ] = DIFFCH;
3958  pHeffect [ CH ] = 0.0;
3959  /* CaCO3 is only mildly soluble */
3960  soluble [ CACO3 ] = 1;
3961  disprob [ CACO3 ] = disbase [ CACO3 ] = 0.10 / DISBIAS;
3962  creates [ CACO3 ] = DIFFCACO3;
3963  pHeffect [ CACO3 ] = 0.0;
3964  /* Slag is not truly soluble, but use its dissolution probability for reaction probability */
3965  soluble [ SLAG ] = 0;
3966  disprob [ SLAG ] = disbase [ SLAG ] = 0.005 / DISBIAS;
3967  creates [ SLAG ] = 0;
3968  pHeffect [ SLAG ] = 1.0;
3969  soluble [ C3AH6 ] = 1;
3970  disprob [ C3AH6 ] = disbase [ C3AH6 ] = 0.01 / DISBIAS; /* changed from 0.5 to 0.01 06.09.00 */
3971  creates [ C3AH6 ] = POROSITY;
3972  pHeffect [ C3AH6 ] = 0.0;
3973  /* Ettringite is initially insoluble */
3974  soluble [ ETTR ] = 0;
3975  /* Changed to 0.008 from 0.020 3/11/99 */
3976  disprob [ ETTR ] = disbase [ ETTR ] = 0.008 / DISBIAS;
3977  creates [ ETTR ] = DIFFETTR;
3978  pHeffect [ ETTR ] = 0.0;
3979  /* Iron-rich ettringite is always insoluble */
3980  soluble [ ETTRC4AF ] = 0;
3981  disprob [ ETTRC4AF ] = disbase [ ETTRC4AF ] = 0.0;
3982  creates [ ETTRC4AF ] = ETTRC4AF;
3983  pHeffect [ ETTRC4AF ] = 0.0;
3984  /* calcium chloride is soluble */
3985  soluble [ CACL2 ] = 1;
3986  disprob [ CACL2 ] = disbase [ CACL2 ] = 0.1 / DISBIAS;
3987  creates [ CACL2 ] = DIFFCACL2;
3988  pHeffect [ CACL2 ] = 0.0;
3989  /* aluminosilicate glass is soluble */
3990  soluble [ ASG ] = 1;
3991  disprob [ ASG ] = disbase [ ASG ] = 0.2 / DISBIAS;
3992  creates [ ASG ] = DIFFAS;
3993  pHeffect [ ASG ] = 1.0;
3994  /* calcium aluminodisilicate is soluble */
3995  soluble [ CAS2 ] = 1;
3996  disprob [ CAS2 ] = disbase [ CAS2 ] = 0.2 / DISBIAS;
3997  creates [ CAS2 ] = DIFFCAS2;
3998  pHeffect [ CAS2 ] = 1.0;
3999 
4000  /* establish molar volumes and heats of formation */
4001  /* molar volumes are in cm^3/mole */
4002  /* heats of formation are in kJ/mole */
4003  /* See paper by Fukuhara et al., Cem. & Conc. Res., 11, 407-14, 1981. */
4004  /* values for Porosity are those of water */
4005  molarv [ POROSITY ] = 18.068;
4006  heatf [ POROSITY ] = ( -285.83 );
4007  waterc [ POROSITY ] = 1.0;
4008  specgrav [ POROSITY ] = 0.99707;
4009 
4010  molarv [ C3S ] = 71.129;
4011  heatf [ C3S ] = ( -2927.82 );
4012  waterc [ C3S ] = 0.0;
4013  specgrav [ C3S ] = 3.21;
4014  /* For improvement in chemical shrinkage correspondence */
4015  /* Changed molar volume of C(1.7)-S-H(4.0) to 108 5/24/95 */
4016  molarv [ CSH ] = 108.;
4017  heatf [ CSH ] = ( -3283.0 );
4018  waterc [ CSH ] = 4.0;
4019  specgrav [ CSH ] = 2.11;
4020 
4021  molarv [ CH ] = 33.1;
4022  heatf [ CH ] = ( -986.1 );
4023  waterc [ CH ] = 1.0;
4024  specgrav [ CH ] = 2.24;
4025 
4026  //uncomment to have the same values as in Pignat's thesis
4027  /*
4028  * molarv[C3S]=72.381;//changed 18.5.05
4029  * specgrav[C3S]=3.15;//changed 18.5.05
4030  * molarv[CSH]=113.5;//changed 18.5.05
4031  * specgrav[CSH]=2.0;//changed 18.5.05
4032  * molarv[CH]=33.036;//changed 18.5.05
4033  */
4034 
4035  /* Assume that calcium carbonate is in the calcite form */
4036  molarv [ CACO3 ] = 36.93;
4037  waterc [ CACO3 ] = 0.0;
4038  specgrav [ CACO3 ] = 2.71;
4039  heatf [ CACO3 ] = ( -1206.92 );
4040 
4041  molarv [ AFMC ] = 261.91;
4042  waterc [ AFMC ] = 11.0;
4043  specgrav [ AFMC ] = 2.17;
4044  /* Need to fill in heat of formation at a later date */
4045  heatf [ AFMC ] = ( 0.0 );
4046 
4047  molarv [ C2S ] = 52.513;
4048  heatf [ C2S ] = ( -2311.6 );
4049  waterc [ C2S ] = 0.0;
4050  specgrav [ C2S ] = 3.28;
4051 
4052  molarv [ C3A ] = 88.94;
4053  heatf [ C3A ] = ( -3587.8 );
4054  waterc [ C3A ] = 0.0;
4055  specgrav [ C3A ] = 3.038;
4056 
4057  molarv [ GYPSUM ] = 74.21;
4058  heatf [ GYPSUM ] = ( -2022.6 );
4059  waterc [ GYPSUM ] = 0.0;
4060  specgrav [ GYPSUM ] = 2.32;
4061  molarv [ GYPSUMS ] = 74.21;
4062  heatf [ GYPSUMS ] = ( -2022.6 );
4063  waterc [ GYPSUMS ] = 0.0;
4064  specgrav [ GYPSUMS ] = 2.32;
4065 
4066  molarv [ ANHYDRITE ] = 52.16;
4067  heatf [ ANHYDRITE ] = ( -1424.6 );
4068  waterc [ ANHYDRITE ] = 0.0;
4069  specgrav [ ANHYDRITE ] = 2.61;
4070 
4071  molarv [ HEMIHYD ] = 52.973;
4072  heatf [ HEMIHYD ] = ( -1574.65 );
4073  waterc [ HEMIHYD ] = 0.0;
4074  specgrav [ HEMIHYD ] = 2.74;
4075 
4076  molarv [ C4AF ] = 130.29;
4077  heatf [ C4AF ] = ( -5090.3 );
4078  waterc [ C4AF ] = 0.0;
4079  specgrav [ C4AF ] = 3.73;
4080 
4081  molarv [ C3AH6 ] = 150.12;
4082  heatf [ C3AH6 ] = ( -5548. );
4083  waterc [ C3AH6 ] = 6.0;
4084  specgrav [ C3AH6 ] = 2.52;
4085 
4086  /* Changed molar volume of FH3 to 69.8 (specific gravity of 3.06) 5/23/95 */
4087  molarv [ FH3 ] = 69.803;
4088  heatf [ FH3 ] = ( -823.9 );
4089  waterc [ FH3 ] = 3.0;
4090  specgrav [ FH3 ] = 3.062;
4091 
4092  molarv [ ETTRC4AF ] = 735.01;
4093  heatf [ ETTRC4AF ] = ( -17539.0 );
4094  waterc [ ETTRC4AF ] = 26.0;
4095  specgrav [ ETTRC4AF ] = 1.7076;
4096  /* Changed molar volue of ettringite to 735 (spec. gr.=1.7076) 5/24/95 */
4097  molarv [ ETTR ] = 735.01;
4098  heatf [ ETTR ] = ( -17539.0 );
4099  waterc [ ETTR ] = 26.0;
4100  specgrav [ ETTR ] = 1.7076;
4101 
4102  molarv [ AFM ] = 312.82;
4103  heatf [ AFM ] = ( -8778.0 );
4104  /* Each mole of AFM which forms requires 12 moles of water, */
4105  /* two of which are supplied by gypsum in forming ETTR */
4106  /* leaving 10 moles to be incorporated from free water */
4107  waterc [ AFM ] = 10.0;
4108  specgrav [ AFM ] = 1.99;
4109 
4110  molarv [ CACL2 ] = 51.62;
4111  heatf [ CACL2 ] = ( -795.8 );
4112  waterc [ CACL2 ] = 0;
4113  specgrav [ CACL2 ] = 2.15;
4114 
4115  molarv [ FREIDEL ] = 296.662;
4116  /* No data available for heat of formation */
4117  heatf [ FREIDEL ] = ( 0.0 );
4118  /* 10 moles of H2O per mole of Freidel's salt */
4119  waterc [ FREIDEL ] = 10.0;
4120  specgrav [ FREIDEL ] = 1.892;
4121 
4122  /* Basic reaction is 2CH + ASG + 6H --> C2ASH8 (Stratlingite) */
4123  molarv [ ASG ] = 49.9;
4124  /* No data available for heat of formation */
4125  heatf [ ASG ] = 0.0;
4126  waterc [ ASG ] = 0.0;
4127  specgrav [ ASG ] = 3.247;
4128 
4129  molarv [ CAS2 ] = 100.62;
4130  /* No data available for heat of formation */
4131  heatf [ CAS2 ] = 0.0;
4132  waterc [ CAS2 ] = 0.0;
4133  specgrav [ CAS2 ] = 2.77;
4134 
4135  molarv [ STRAT ] = 215.63;
4136  /* No data available for heat of formation */
4137  heatf [ STRAT ] = 0.0;
4138  /* 8 moles of water per mole of stratlingite */
4139  waterc [ STRAT ] = 8.0;
4140  specgrav [ STRAT ] = 1.94;
4141 
4142  molarv [ POZZ ] = 27.0;
4143  /* Use heat of formation of SiO2 (quartz) for unreacted pozzolan */
4144  /* Source- Handbook of Chemistry and Physics */
4145  heatf [ POZZ ] = -907.5;
4146  waterc [ POZZ ] = 0.0;
4147  specgrav [ POZZ ] = 2.22;
4148 
4149  /* Data for Pozzolanic CSH based on work of Atlassi, DeLarrard, */
4150  /* and Jensen */
4151  /* gives a chemical shrinkage of 0.2 g H2O/g CSF */
4152  /* heat of formation estimated based on heat release of */
4153  /* 780 J/g Condensed Silica Fume */
4154  /* Changed stoichiometry to be C(1.1)SH(3.9) to see effect on */
4155  /* results 1/22/97 */
4156  /* MW is 191.8 g/mole */
4157  /* Changed molar volume to 101.81 3/10/97 */
4158  molarv [ POZZCSH ] = 101.81;
4159  waterc [ POZZCSH ] = 3.9;
4160  specgrav [ POZZCSH ] = 1.884;
4161  heatf [ POZZCSH ] = ( -2.1 );
4162 
4163  /* Assume inert filler has same specific gravity and molar volume as SiO2 */
4164  molarv [ INERT ] = 27.0;
4165  heatf [ INERT ] = 0.0;
4166  waterc [ INERT ] = 0.0;
4167  specgrav [ INERT ] = 2.2;
4168 
4169  molarv [ ABSGYP ] = 74.21;
4170  heatf [ ABSGYP ] = ( -2022.6 );
4171  waterc [ ABSGYP ] = 0.0;
4172  specgrav [ ABSGYP ] = 2.32;
4173 
4174  molarv [ EMPTYP ] = 18.068;
4175  heatf [ EMPTYP ] = ( -285.83 );
4176  waterc [ EMPTYP ] = 0.0;
4177  specgrav [ EMPTYP ] = 0.99707;
4178 
4179  /* Read in values for alkali characteristics and */
4180  /* convert them to fractions from percentages */
4181  // First define non-zero values due to divisions
4182  //totsodium=0.2/100;
4183  //totpotassium=0.5/100;
4184  //rssodium=0.03/100;
4185  //rspotassium=0.25/100;
4186  specgrav [ SLAG ] = 2.87;
4187  specgrav [ SLAGCSH ] = 2.35;
4188  molarv [ SLAG ] = 945.72;
4189  molarv [ SLAGCSH ] = 1717.53;
4190  slagcasi = 1.433;
4191  slaghydcasi = 1.35;
4192  siperslag = 15.0;
4193  waterc [ SLAGCSH ] = 5.533 * siperslag;
4194  slagc3a = 1.0;
4195  slagreact = 1.0;
4196 
4197 #ifdef CMLFILE
4198  F->get_value(19, totsodium);
4199  F->get_value(20, totpotassium);
4200  F->get_value(21, rssodium);
4201  F->get_value(22, rspotassium);
4202 #endif
4203 #ifdef TINYXML
4204  QueryNumAttributeExt(xmlFile, "Total_sodium", 0, totsodium);
4205  QueryNumAttributeExt(xmlFile, "Total_potassium", 0, totpotassium);
4206  QueryNumAttributeExt(xmlFile, "Readily_soluble_sodium", 0, rssodium);
4207  QueryNumAttributeExt(xmlFile, "Readily_soluble_potassium", 0, rspotassium);
4208 #endif
4209 
4210  /*fscanf(alkalifile,"%f",&totsodium);
4211  * fscanf(alkalifile,"%f",&totpotassium);
4212  * fscanf(alkalifile,"%f",&rssodium);
4213  * fscanf(alkalifile,"%f",&rspotassium);
4214  */
4215 
4216  totsodium /= 100.;
4217  totpotassium /= 100.;
4218  rssodium /= 100.;
4219  rspotassium /= 100.;
4220 
4221 
4222  if ( pHactive ) {
4223  /* Read in values for slag characteristics */
4224 
4225  if ( ( slagfile = fopen("slagchar.dat", "r") ) == NULL ) {
4226  OOFEM_ERROR("Slag file can not be opened");
4227  }
4228 
4229 
4230  if ( fscanf(slagfile, "%lf", & slagin) != 1 ) {
4231  OOFEM_ERROR("slagfile reading error");
4232  }
4233 
4234  if ( fscanf(slagfile, "%lf", & slagin) != 1 ) {
4235  OOFEM_ERROR("slagfile reading error");
4236  }
4237 
4238  if ( fscanf(slagfile, "%lf", & slagin) != 1 ) {
4239  OOFEM_ERROR("slagfile reading error");
4240  }
4241 
4242  specgrav [ SLAG ] = slagin;
4243  if ( fscanf(slagfile, "%lf", & slagin) != 1 ) {
4244  OOFEM_ERROR("slagfile reading error");
4245  }
4246 
4247  specgrav [ SLAGCSH ] = slagin;
4248  if ( fscanf(slagfile, "%lf", & slagin) != 1 ) {
4249  OOFEM_ERROR("slagfile reading error");
4250  }
4251 
4252  molarv [ SLAG ] = slagin;
4253  if ( fscanf(slagfile, "%lf", & slagin) != 1 ) {
4254  OOFEM_ERROR("slagfile reading error");
4255  }
4256 
4257  molarv [ SLAGCSH ] = slagin;
4258  if ( fscanf(slagfile, "%lf", & slagcasi) != 1 ) {
4259  OOFEM_ERROR("slagfile reading error");
4260  }
4261 
4262  if ( fscanf(slagfile, "%lf", & slaghydcasi) != 1 ) {
4263  OOFEM_ERROR("slagfile reading error");
4264  }
4265 
4266  if ( fscanf(slagfile, "%lf", & siperslag) != 1 ) {
4267  OOFEM_ERROR("slagfile reading error");
4268  }
4269 
4270  if ( fscanf(slagfile, "%lf", & slagin) != 1 ) {
4271  OOFEM_ERROR("slagfile reading error");
4272  }
4273 
4274  waterc [ SLAGCSH ] = slagin * siperslag;
4275  if ( fscanf(slagfile, "%lf", & slagc3a) != 1 ) {
4276  OOFEM_ERROR("slagfile reading error");
4277  }
4278 
4279  if ( fscanf(slagfile, "%lf", & slagreact) != 1 ) {
4280  OOFEM_ERROR("slagfile reading error");
4281  }
4282 
4283  waterc [ SLAG ] = 0.0;
4284  heatf [ SLAG ] = 0.0;
4285  heatf [ SLAGCSH ] = 0.0;
4286  fclose(slagfile);
4287  }
4288 
4289  //use default values
4290  specgrav [ SLAG ] = 2.87;
4291  specgrav [ SLAGCSH ] = 2.35;
4292  molarv [ SLAG ] = 945.72;
4293  molarv [ SLAGCSH ] = 1717.53;
4294  slagcasi = 1.4333;
4295  slaghydcasi = 1.35;
4296  siperslag = 15.0;
4297  waterc [ SLAGCSH ] = 5.533 * siperslag;
4298  slagc3a = 1.0;
4299  slagreact = 1.0;
4300 
4301  /* Compute slag probabilities as defined above */
4302  CHperslag = siperslag * ( slaghydcasi - slagcasi ) + 3. * slagc3a;
4303  if ( CHperslag < 0.0 ) {
4304  CHperslag = 0.0;
4305  }
4306 
4307  p2slag = ( molarv [ SLAG ] + molarv [ CH ] * CHperslag + molarv [ POROSITY ] * ( waterc [ SLAGCSH ] - CHperslag + waterc [ C3AH6 ] * slagc3a ) - molarv [ SLAGCSH ] - molarv [ C3AH6 ] * slagc3a ) / molarv [ SLAG ];
4308  p1slag = 1.0 - p2slag;
4309  p3slag = ( molarv [ SLAGCSH ] / molarv [ SLAG ] ) - p1slag;
4310  p4slag = CHperslag * molarv [ CH ] / molarv [ SLAG ];
4311  p5slag = slagc3a * molarv [ C3A ] / molarv [ SLAG ];
4312  if ( p5slag > 1.0 ) {
4313  p5slag = 1.0;
4314  printf("Error in range of C3A/slag value... reset to 1.0 \n");
4315  }
4316 }
4317 
4318 /* routine to check if a pixel located at (xck,yck,zck) is on an edge */
4319 /* (in contact with pore space) in 3-D system */
4320 /* Called by passone */
4321 /* Calls no other routines */
4322 int CemhydMatStatus :: chckedge(int xck, int yck, int zck)
4323 {
4324  int edgeback, x2, y2, z2;
4325  int ip;
4326 
4327  edgeback = 0;
4328 
4329  /* Check all six neighboring pixels */
4330  /* with periodic boundary conditions */
4331  for ( ip = 0; ( ( ip < NEIGHBORS ) && ( edgeback == 0 ) ); ip++ ) {
4332  x2 = xck + xoff [ ip ];
4333  y2 = yck + yoff [ ip ];
4334  z2 = zck + zoff [ ip ];
4335  if ( x2 >= SYSIZE ) {
4336  x2 = 0;
4337  }
4338 
4339  if ( x2 < 0 ) {
4340  x2 = SYSIZEM1;
4341  }
4342 
4343  if ( y2 >= SYSIZE ) {
4344  y2 = 0;
4345  }
4346 
4347  if ( y2 < 0 ) {
4348  y2 = SYSIZEM1;
4349  }
4350 
4351  if ( z2 >= SYSIZE ) {
4352  z2 = 0;
4353  }
4354 
4355  if ( z2 < 0 ) {
4356  z2 = SYSIZEM1;
4357  }
4358 
4359  if ( mic [ x2 ] [ y2 ] [ z2 ] == POROSITY ) {
4360  edgeback = 1;
4361  }
4362  }
4363 
4364  return ( edgeback );
4365 }
4366 
4367 /* routine for first pass through microstructure during dissolution */
4368 /* low and high indicate phase ID range to check for surface sites */
4369 /* Called by dissolve */
4370 /* Calls chckedge */
4371 void CemhydMatStatus :: passone(int low, int high, int cycid, int cshexflag)
4372 {
4373  int i, xid, yid, zid, phid, edgef, phread, cshcyc;
4374 
4375  /* gypready used to determine if any soluble gypsum remains */
4376  if ( ( low <= GYPSUM ) && ( GYPSUM <= high ) ) {
4377  gypready = 0;
4378  }
4379 
4380  /* Zero out count for the relevant phases */
4381  for ( i = low; i <= high; i++ ) {
4382  count [ i ] = 0;
4383  }
4384 
4385  /* Scan the entire 3-D microstructure */
4386  for ( xid = 0; xid < SYSIZE; xid++ ) {
4387  for ( yid = 0; yid < SYSIZE; yid++ ) {
4388  for ( zid = 0; zid < SYSIZE; zid++ ) {
4389  phread = mic [ xid ] [ yid ] [ zid ];
4390  /* Update heat data and water consumed for solid CSH */
4391  if ( ( cshexflag == 1 ) && ( phread == CSH ) ) {
4392  cshcyc = cshage [ xid ] [ yid ] [ zid ];
4393  heatsum += heatf [ CSH ] / molarvcsh [ cshcyc ];
4394  molesh2o += watercsh [ cshcyc ] / molarvcsh [ cshcyc ];
4395  }
4396 
4397  /* Identify phase and update count */
4398  phid = 60;
4399  for ( i = low; ( ( i <= high ) && ( phid == 60 ) ); i++ ) {
4400  if ( mic [ xid ] [ yid ] [ zid ] == i ) {
4401  phid = i;
4402  /* Update count for this phase */
4403  count [ i ] += 1;
4404  if ( ( i == GYPSUM ) || ( i == GYPSUMS ) ) {
4405  gypready += 1;
4406  }
4407 
4408  /* If first cycle, then accumulate initial counts */
4409  if ( cycid == 1 ) { //fixed (ncyc cancelled)
4410  if ( i == POROSITY ) {
4411  porinit += 1;
4412  }
4413  /* Ordered in terms of likely volume fractions */
4414  /* (largest to smallest) to speed execution */
4415  else if ( i == C3S ) {
4416  c3sinit += 1;
4417  } else if ( i == C2S ) {
4418  c2sinit += 1;
4419  } else if ( i == C3A ) {
4420  c3ainit += 1;
4421  } else if ( i == C4AF ) {
4422  c4afinit += 1;
4423  } else if ( i == GYPSUM ) {
4424  ncsbar += 1;
4425  } else if ( i == GYPSUMS ) {
4426  ncsbar += 1;
4427  } else if ( i == ANHYDRITE ) {
4428  anhinit += 1;
4429  } else if ( i == HEMIHYD ) {
4430  heminit += 1;
4431  } else if ( i == POZZ ) {
4432  nfill += 1;
4433  } else if ( i == SLAG ) {
4434  slaginit += 1;
4435  } else if ( i == ETTR ) {
4436  netbar += 1;
4437  } else if ( i == ETTRC4AF ) {
4438  netbar += 1;
4439  }
4440  }
4441  }
4442  }
4443 
4444  if ( phid != 60 ) {
4445  /* If phase is soluble, see if it is in contact with porosity */
4446  if ( ( cycid != 0 ) && ( soluble [ phid ] == 1 ) ) {
4447  edgef = chckedge(xid, yid, zid);
4448  if ( edgef == 1 ) {
4449  /* Surface eligible species has an ID OFFSET greater than its original value */
4450  mic [ xid ] [ yid ] [ zid ] += OFFSET;
4451  }
4452  }
4453  }
4454  }
4455 
4456  /* end of zid */
4457  }
4458 
4459  /* end of yid */
4460  }
4461 
4462  /* end of xid */
4463 }
4464 
4465 /* routine to locate a diffusing CSH species near dissolution source */
4466 /* at (xcur,ycur,zcur) */
4467 /* Called by dissolve */
4468 /* Calls no other routines */
4469 int CemhydMatStatus :: loccsh(int xcur, int ycur, int zcur, int extent)
4470 {
4471  int effort, tries, xmod, ymod, zmod;
4472  struct ants *antnew;
4473 
4474  effort = 0; /* effort indicates if appropriate location found */
4475  tries = 0;
4476  /* 500 tries in immediate vicinity */
4477  while ( ( effort == 0 ) && ( tries < 500 ) ) {
4478  tries += 1;
4479  xmod = ( -extent ) + ( int ) ( ( 2. * ( float ) extent + 1. ) * ran1(seed) );
4480  ymod = ( -extent ) + ( int ) ( ( 2. * ( float ) extent + 1. ) * ran1(seed) );
4481  zmod = ( -extent ) + ( int ) ( ( 2. * ( float ) extent + 1. ) * ran1(seed) );
4482  if ( xmod > extent ) {
4483  xmod = extent;
4484  }
4485 
4486  if ( ymod > extent ) {
4487  ymod = extent;
4488  }
4489 
4490  if ( zmod > extent ) {
4491  zmod = extent;
4492  }
4493 
4494  xmod += xcur;
4495  ymod += ycur;
4496  zmod += zcur;
4497  /* Periodic boundaries */
4498  if ( xmod >= SYSIZE ) {
4499  xmod -= SYSIZE;
4500  } else if ( xmod < 0 ) {
4501  xmod += SYSIZE;
4502  }
4503 
4504  if ( zmod >= SYSIZE ) {
4505  zmod -= SYSIZE;
4506  } else if ( zmod < 0 ) {
4507  zmod += SYSIZE;
4508  }
4509 
4510  if ( ymod < 0 ) {
4511  ymod += SYSIZE;
4512  } else if ( ymod >= SYSIZE ) {
4513  ymod -= SYSIZE;
4514  }
4515 
4516  if ( mic [ xmod ] [ ymod ] [ zmod ] == POROSITY ) {
4517  effort = 1;
4518  mic [ xmod ] [ ymod ] [ zmod ] = DIFFCSH;
4519  nmade += 1;
4520  ngoing += 1;
4521  /* Add CSH diffusing species to the linked list */
4522  antnew = ( struct ants * ) malloc( sizeof( struct ants ) );
4523  antnew->x = xmod;
4524  antnew->y = ymod;
4525  antnew->z = zmod;
4526  antnew->id = DIFFCSH;
4527  antnew->cycbirth = cyccnt;
4528  /* Now connect this ant structure to end of linked list */
4529  antnew->prevant = tailant;
4530  tailant->nextant = antnew;
4531  antnew->nextant = NULL;
4532  tailant = antnew;
4533  }
4534  }
4535 
4536  return ( effort );
4537 }
4538 
4539 /* routine to count number of pore pixels in a cube of size boxsize */
4540 /* centered at (qx,qy,qz) */
4541 /* Called by makeinert */
4542 /* Calls no other routines */
4543 int CemhydMatStatus :: countbox(int boxsize, int qx, int qy, int qz)
4544 {
4545  int nfound, ix, iy, iz, qxlo, qxhi, qylo, qyhi, qzlo, qzhi;
4546  int hx, hy, hz, boxhalf;
4547 
4548  boxhalf = boxsize / 2;
4549  nfound = 0;
4550  qxlo = qx - boxhalf;
4551  qxhi = qx + boxhalf;
4552  qylo = qy - boxhalf;
4553  qyhi = qy + boxhalf;
4554  qzlo = qz - boxhalf;
4555  qzhi = qz + boxhalf;
4556  /* Count the number of requisite pixels in the 3-D cube box */
4557  /* using periodic boundaries */
4558  for ( ix = qxlo; ix <= qxhi; ix++ ) {
4559  hx = ix;
4560  if ( hx < 0 ) {
4561  hx += SYSIZE;
4562  } else if ( hx >= SYSIZE ) {
4563  hx -= SYSIZE;
4564  }
4565 
4566  for ( iy = qylo; iy <= qyhi; iy++ ) {
4567  hy = iy;
4568  if ( hy < 0 ) {
4569  hy += SYSIZE;
4570  } else if ( hy >= SYSIZE ) {
4571  hy -= SYSIZE;
4572  }
4573 
4574  for ( iz = qzlo; iz <= qzhi; iz++ ) {
4575  hz = iz;
4576  if ( hz < 0 ) {
4577  hz += SYSIZE;
4578  } else if ( hz >= SYSIZE ) {
4579  hz -= SYSIZE;
4580  }
4581 
4582  /* Count if porosity, diffusing species, or empty porosity */
4583  if ( ( mic [ hx ] [ hy ] [ hz ] < C3S ) || ( mic [ hx ] [ hy ] [ hz ] > ABSGYP ) ) {
4584  nfound += 1;
4585  }
4586  }
4587  }
4588  }
4589 
4590  return ( nfound );
4591 }
4592 
4593 /* routine to count number of special pixels in a cube of size boxsize */
4594 /* centered at (qx,qy,qz) */
4595 /* special pixels are those not belonging to one of the cement clinker, */
4596 /* calcium sulfate, or pozzolanic mineral admixture phases */
4597 /* Called by addrand */
4598 /* Calls no other routines */
4599 int CemhydMatStatus :: countboxc(int boxsize, int qx, int qy, int qz)
4600 {
4601  int nfound, ix, iy, iz, qxlo, qxhi, qylo, qyhi, qzlo, qzhi;
4602  int hx, hy, hz, boxhalf;
4603 
4604  boxhalf = boxsize / 2;
4605  nfound = 0;
4606  qxlo = qx - boxhalf;
4607  qxhi = qx + boxhalf;
4608  qylo = qy - boxhalf;
4609  qyhi = qy + boxhalf;
4610  qzlo = qz - boxhalf;
4611  qzhi = qz + boxhalf;
4612  /* Count the number of requisite pixels in the 3-D cube box */
4613  /* using periodic boundaries */
4614  for ( ix = qxlo; ix <= qxhi; ix++ ) {
4615  hx = ix;
4616  if ( hx < 0 ) {
4617  hx += SYSIZE;
4618  } else if ( hx >= SYSIZE ) {
4619  hx -= SYSIZE;
4620  }
4621 
4622  for ( iy = qylo; iy <= qyhi; iy++ ) {
4623  hy = iy;
4624  if ( hy < 0 ) {
4625  hy += SYSIZE;
4626  } else if ( hy >= SYSIZE ) {
4627  hy -= SYSIZE;
4628  }
4629 
4630  for ( iz = qzlo; iz <= qzhi; iz++ ) {
4631  hz = iz;
4632  if ( hz < 0 ) {
4633  hz += SYSIZE;
4634  } else if ( hz >= SYSIZE ) {
4635  hz -= SYSIZE;
4636  }
4637 
4638  /* Count if not cement clinker */
4639  if ( ( mic [ hx ] [ hy ] [ hz ] < C3S ) || ( mic [ hx ] [ hy ] [ hz ] > POZZ ) ) {
4640  nfound += 1;
4641  }
4642  }
4643  }
4644  }
4645 
4646  return ( nfound );
4647 }
4648 
4649 /* routine to create ndesire pixels of empty pore space to simulate */
4650 /* self-desiccation, or external drying */
4651 /* Called by dissolve */
4652 /* Calls countbox */
4653 /* each ....togo structure contains information only at x,y,z*/
4654 /* headtogo is the first voxel where the water is taken from*/
4655 
4656 
4657 void CemhydMatStatus :: makeinert(long int ndesire) {
4658  struct togo *headtogo, *tailtogo, *newtogo, *lasttogo, *onetogo;
4659  long int idesire;
4660  int px, py, pz, placed, cntpore, cntmax;
4661 
4662  /* First allocate the first element of the linked list */
4663  headtogo = ( struct togo * ) malloc( sizeof( struct togo ) );
4664  headtogo->x = headtogo->y = headtogo->z = ( -1 );
4665  headtogo->npore = 0;
4666  headtogo->nexttogo = NULL;
4667  headtogo->prevtogo = NULL;
4668  tailtogo = headtogo;
4669  cntmax = 0;
4670 
4671 #ifdef PRINTF
4672  printf("In makeinert with %ld needed elements \n", ndesire);
4673 #endif
4674 
4675  fflush(stdout);
4676  /* Add needed number of elements to the end of the list */
4677  for ( idesire = 2; idesire <= ndesire; idesire++ ) {
4678  newtogo = ( struct togo * ) malloc( sizeof( struct togo ) );
4679  newtogo->npore = 0;
4680  newtogo->x = newtogo->y = newtogo->z = ( -1 );
4681  tailtogo->nexttogo = newtogo;
4682  newtogo->prevtogo = tailtogo;
4683  tailtogo = newtogo;
4684  }
4685 
4686  /* Now scan the microstructure and rank the sites */
4687  for ( px = 0; px < SYSIZE; px++ ) {
4688  for ( py = 0; py < SYSIZE; py++ ) {
4689  for ( pz = 0; pz < SYSIZE; pz++ ) {
4690  if ( mic [ px ] [ py ] [ pz ] == POROSITY ) {
4691  cntpore = countbox(cubesize, px, py, pz);
4692  if ( cntpore > cntmax ) {
4693  cntmax = cntpore;
4694  }
4695 
4696  /* Store this site value at appropriate place in */
4697  /* sorted linked list */
4698  if ( cntpore > ( tailtogo->npore ) ) {
4699  placed = 0;
4700  lasttogo = tailtogo;
4701  while ( placed == 0 ) {
4702  newtogo = lasttogo->prevtogo;
4703  if ( newtogo == NULL ) { //if at the beginning of the list (tailtogo->prevtogo==NULL)
4704  placed = 2;
4705  } else {
4706  if ( cntpore <= ( newtogo->npore ) ) {
4707  placed = 1;
4708  }
4709  }
4710 
4711  if ( placed == 0 ) {
4712  lasttogo = newtogo;
4713  }
4714  }
4715 
4716  onetogo = ( struct togo * ) malloc( sizeof( struct togo ) );
4717  onetogo->x = px;
4718  onetogo->y = py;
4719  onetogo->z = pz;
4720  onetogo->npore = cntpore;
4721  /* Insertion at the head of the list */
4722  if ( placed == 2 ) {
4723  onetogo->prevtogo = NULL;
4724  onetogo->nexttogo = headtogo;
4725  headtogo->prevtogo = onetogo;
4726  headtogo = onetogo;
4727  }
4728 
4729  if ( placed == 1 ) {
4730  onetogo->nexttogo = lasttogo;
4731  onetogo->prevtogo = newtogo;
4732  lasttogo->prevtogo = onetogo;
4733  newtogo->nexttogo = onetogo;
4734  }
4735 
4736  /* Eliminate the last element */
4737  lasttogo = tailtogo;
4738  tailtogo = tailtogo->prevtogo;
4739  tailtogo->nexttogo = NULL;
4740  //printf("DEMEM lasttogo1 addr %p\n", lasttogo);
4741  free(lasttogo);
4742  }
4743  }
4744  }
4745  }
4746  }
4747 
4748  /* Now remove the sites */
4749  /* starting at the head of the list */
4750  /* and deallocate all of the used memory */
4751  for ( idesire = 1; idesire <= ndesire; idesire++ ) {
4752  px = headtogo->x;
4753  py = headtogo->y;
4754  pz = headtogo->z;
4755  if ( px != ( -1 ) ) {
4756  mic [ px ] [ py ] [ pz ] = EMPTYP;
4757  count [ POROSITY ] -= 1;
4758  count [ EMPTYP ] += 1;
4759  }
4760 
4761  lasttogo = headtogo;
4762  headtogo = headtogo->nexttogo;
4763  //printf("DEMEM lasttogo2 addr %p\n", (const void *)lasttogo);
4764  free(lasttogo);
4765  }
4766 
4767  /* If only small cubes of porosity were found, then adjust */
4768  /* cubesize to have a more efficient search in the future */
4769  if ( cubesize > CUBEMIN ) {
4770  if ( ( 2 * cntmax ) < ( cubesize * cubesize * cubesize ) ) {
4771  cubesize -= 2;
4772  }
4773  }
4774 }
4775 
4776 
4777 /* routine to add extra SLAG CSH when SLAG reacts */
4778 /* SLAG located at (xpres,ypres,zpres) */
4779 /* Called by dissolve */
4780 /* Calls moveone and edgecnt */
4781 void CemhydMatStatus :: extslagcsh(int xpres, int ypres, int zpres) {
4782  int check, sump, xchr, ychr, zchr, fchr, i1, action, numnear;
4783  long int tries;
4784  int mstest = 0, mstest2 = 0;
4785 
4786  /* first try 6 neighboring locations until */
4787  /* a) successful */
4788  /* b) all 6 sites are tried or */
4789  /* c) 100 tries are made */
4790  /* try to grow slag C-S-H as plates */
4791  fchr = 0;
4792  sump = 1;
4793  for ( i1 = 1; ( ( i1 <= 100 ) && ( fchr == 0 ) && ( sump != 30030 ) ); i1++ ) {
4794  /* determine location of neighbor (using periodic boundaries) */
4795  xchr = xpres;
4796  ychr = ypres;
4797  zchr = zpres;
4798  action = 0;
4799  sump *= moveone(& xchr, & ychr, & zchr, & action, sump);
4800  if ( action == 0 ) {
4801  printf("Error in value of action in extpozz \n");
4802  }
4803 
4804  check = mic [ xchr ] [ ychr ] [ zchr ];
4805  /* Determine the direction of the neighbor selected and */
4806  /* the plates possible for growth */
4807  if ( xchr != xpres ) {
4808  mstest = 1;
4809  mstest2 = 2;
4810  }
4811 
4812  if ( ychr != ypres ) {
4813  mstest = 2;
4814  mstest2 = 3;
4815  }
4816 
4817  if ( zchr != zpres ) {
4818  mstest = 3;
4819  mstest2 = 1;
4820  }
4821 
4822  /* if neighbor is porosity, locate the SLAG CSH there */
4823  if ( check == POROSITY ) {
4824  if ( ( faces [ xpres ] [ ypres ] [ zpres ] == 0 ) || ( mstest == faces [ xpres ] [ ypres ] [ zpres ] ) || ( mstest2 == faces [ xpres ] [ ypres ] [ zpres ] ) ) {
4825  mic [ xchr ] [ ychr ] [ zchr ] = SLAGCSH;
4826  faces [ xchr ] [ ychr ] [ zchr ] = faces [ xpres ] [ ypres ] [ zpres ];
4827  count [ SLAGCSH ] += 1;
4828  count [ POROSITY ] -= 1;
4829  fchr = 1;
4830  }
4831  }
4832  }
4833 
4834  /* if no neighbor available, locate SLAG CSH at random location */
4835  /* in pore space */
4836  tries = 0;
4837  while ( fchr == 0 ) {
4838  tries += 1;
4839  /* generate a random location in the 3-D system */
4840  xchr = ( int ) ( ( float ) SYSIZE * ran1(seed) );
4841  ychr = ( int ) ( ( float ) SYSIZE * ran1(seed) );
4842  zchr = ( int ) ( ( float ) SYSIZE * ran1(seed) );
4843  if ( xchr >= SYSIZE ) {
4844  xchr = 0;
4845  }
4846 
4847  if ( ychr >= SYSIZE ) {
4848  ychr = 0;
4849  }
4850 
4851  if ( zchr >= SYSIZE ) {
4852  zchr = 0;
4853  }
4854 
4855  check = mic [ xchr ] [ ychr ] [ zchr ];
4856  /* if location is porosity, locate the extra SLAG CSH there */
4857  if ( check == POROSITY ) {
4858  numnear = edgecnt(xchr, ychr, zchr, SLAG, CSH, SLAGCSH);
4859  /* Be sure that one neighboring species is CSH or */
4860  /* SLAG material */
4861  if ( ( tries > 5000 ) || ( numnear < 26 ) ) {
4862  mic [ xchr ] [ ychr ] [ zchr ] = SLAGCSH;
4863  count [ SLAGCSH ] += 1;
4864  count [ POROSITY ] -= 1;
4865  fchr = 1;
4866  }
4867  }
4868  }
4869 }
4870 
4871 /* routine to implement a cycle of dissolution */
4872 /* Called by main program */
4873 /* Calls passone, loccsh, and makeinert */
4875 {
4876  int nc3aext, ncshext, nchext, ngypext, nanhext, plok;
4877  int nsum5, nsum4, nsum3, nsum2, nhemext, nsum6, nc4aext;
4878  int phid, phnew, plnew, cread;
4879  int i, xloop, yloop, zloop, xc, yc;
4880  int zc, cycnew;
4881  long int ctest;
4882  long int tries;
4883  int placed, cshrand, maxsulfate, msface;
4884  long int ncshgo, nsurf, suminit;
4885  long int xext, nhgd, npchext, nslagc3a = 0;
4886  float pdis, plfh3, fchext, fc3aext, fanhext, mass_now, mass_fa_now, tot_mass, heatfill;
4887  float dfact, dfact1, molesdh2o, h2oinit, heat4, fhemext, fc4aext;
4888  float pconvert, pc3scsh, pc2scsh, calcx, calcy, calcz, tdisfact;
4889  float frafm, frettr, frhyg, frtot, mc3ar, mc4ar, p3init;
4890  struct ants *antadd;
4891 
4892  /* Initialize variables */
4893  nmade = 0;
4894  npchext = ncshgo = cshrand = 0; /* counter for number of csh diffusing species to */
4895  /* be located at random locations in microstructure */
4896  heat_old = heat_new; /* new and old values for heat released */
4897 
4898  /* Initialize dissolution and phase counters */
4899  nsurf = 0;
4900  for ( i = 0; i <= EMPTYP; i++ ) {
4901  discount [ i ] = 0;
4902  count [ i ] = 0;
4903  }
4904 
4905  /* Pass one- highlight all edge points which are soluble */
4906  soluble [ C3AH6 ] = 0;
4907  heatsum = 0.0;
4908  molesh2o = 0.0;
4909  passone(0, EMPTYP, cycle, 1);
4910 #ifdef PRINTF
4911  printf("Returned from passone \n");
4912 #endif
4913  fflush(stdout);
4914 
4915  sulf_solid = count [ GYPSUM ] + count [ GYPSUMS ] + count [ HEMIHYD ] + count [ ANHYDRITE ];
4916  /* If first cycle, then determine all mixture proportions based */
4917  /* on user input and original microstructure */
4918  if ( cycle == 1 ) {
4919  /* Mass of cement in system */
4920  cemmass = ( specgrav [ C3S ] * ( float ) count [ C3S ] + specgrav [ C2S ] *
4921  ( float ) count [ C2S ] + specgrav [ C3A ] * ( float ) count [ C3A ] +
4922  specgrav [ C4AF ] * ( float ) count [ C4AF ] );
4923  /* +specgrav[GYPSUM]*
4924  * (float)count[GYPSUM]+specgrav[ANHYDRITE]*(float)
4925  * count[ANHYDRITE]+specgrav[HEMIHYD]*(float)count[HEMIHYD]); */
4926  cemmasswgyp = ( specgrav [ C3S ] * ( float ) count [ C3S ] + specgrav [ C2S ] *
4927  ( float ) count [ C2S ] + specgrav [ C3A ] * ( float ) count [ C3A ] +
4928  specgrav [ C4AF ] * ( float ) count [ C4AF ] + specgrav [ GYPSUM ] *
4929  ( float ) count [ GYPSUM ] + specgrav [ ANHYDRITE ] * ( float )
4930  count [ ANHYDRITE ] + specgrav [ HEMIHYD ] * ( float ) count [ HEMIHYD ] );
4931  flyashmass = ( specgrav [ ASG ] * ( float ) count [ ASG ] + specgrav [ CAS2 ] *
4932  ( float ) count [ CAS2 ] + specgrav [ POZZ ] * ( float ) count [ POZZ ] );
4933  CH_mass = specgrav [ CH ] * ( float ) count [ CH ];
4934  /* Total mass in system neglecting single aggregate */
4935  tot_mass = cemmass + ( float ) count [ POROSITY ] + specgrav [ INERT ] *
4936  ( float ) count [ INERT ] + specgrav [ CACL2 ] * ( float ) count [ CACL2 ] +
4937  specgrav [ ASG ] * ( float ) count [ ASG ] +
4938  specgrav [ SLAG ] * ( float ) count [ SLAG ] +
4939  specgrav [ HEMIHYD ] * ( float ) count [ HEMIHYD ] +
4940  specgrav [ ANHYDRITE ] * ( float ) count [ ANHYDRITE ] +
4941  specgrav [ CAS2 ] * ( float ) count [ CAS2 ] +
4942  specgrav [ CSH ] * ( float ) count [ CSH ] +
4943  specgrav [ GYPSUM ] * ( float ) count [ GYPSUM ] +
4944  specgrav [ ANHYDRITE ] * ( float ) count [ ANHYDRITE ] +
4945  specgrav [ HEMIHYD ] * ( float ) count [ HEMIHYD ] +
4946  specgrav [ GYPSUMS ] * ( float ) count [ GYPSUMS ] +
4947  specgrav [ POZZ ] * ( float ) count [ POZZ ] + CH_mass;
4948  /* water-to-cement ratio */
4949  if ( cemmass != 0.0 ) {
4950  w_to_c = ( float ) count [ POROSITY ] / ( cemmass +
4951  specgrav [ GYPSUM ] * ( float ) count [ GYPSUM ] +
4952  specgrav [ ANHYDRITE ] * ( float ) count [ ANHYDRITE ] +
4953  specgrav [ HEMIHYD ] * ( float ) count [ HEMIHYD ] );
4954  } else {
4955  w_to_c = 0.0;
4956  }
4957 
4958  /* totfract is the total cement volume count including calcium sulfates */
4959  /* fractwithfill is the total count of cement and solid fillers and */
4960  /* mineral admixtures DPB- 10/04 */
4961  totfract = count [ C3S ] + count [ C2S ];
4962  totfract += ( count [ C3A ] + count [ C4AF ] );
4963  totfract += ( count [ GYPSUM ] + count [ ANHYDRITE ] + count [ HEMIHYD ] );
4964  fractwithfill = totfract + count [ CACO3 ] + count [ SLAG ] + count [ INERT ];
4965  fractwithfill += count [ POZZ ] + count [ CAS2 ] + count [ ASG ] + count [ CACL2 ];
4966  totfract /= ( float ) SYSIZE;
4967  totfract /= ( float ) SYSIZE;
4968  totfract /= ( float ) SYSIZE;
4969  fractwithfill /= ( float ) SYSIZE;
4970  fractwithfill /= ( float ) SYSIZE;
4971  fractwithfill /= ( float ) SYSIZE;
4972  /* Adjust masses for presence of aggregates in concrete */
4973  mass_water = ( 1. - mass_agg ) * ( float ) count [ POROSITY ] / tot_mass;
4974  mass_CH = ( 1. - mass_agg ) * CH_mass / tot_mass;
4975  /* pozzolan-to-cement ratio */
4976  if ( cemmass != 0.0 ) {
4977  s_to_c = ( float ) ( count [ INERT ] * specgrav [ INERT ] +
4978  count [ CACL2 ] * specgrav [ CACL2 ] + count [ ASG ] * specgrav [ ASG ] +
4979  count [ CAS2 ] * specgrav [ CAS2 ] +
4980  count [ SLAG ] * specgrav [ SLAG ] +
4981  count [ POZZ ] * specgrav [ POZZ ] ) / cemmass;
4982  } else {
4983  s_to_c = 0.0;
4984  }
4985 
4986  /* Conversion factor to kJ/kg for heat produced */
4987  if ( cemmass != 0.0 ) {
4988  heatfill = ( float ) ( count [ INERT ] + count [ SLAG ] + count [ POZZ ] + count [ CACL2 ] + count [ ASG ] + count [ CAS2 ] ) / cemmass; //units pixels^3/g, inverse to the density
4989  } else {
4990  heatfill = 0.0;
4991  }
4992 
4993  if ( w_to_c > 0.01 ) {
4994  heat_cf = 1000. / SYSIZE_POW3 * ( 0.3125 + w_to_c + heatfill ); //heat conversion factor,1/ro_c=0.3125, fixed
4995  } else {
4996  /* Need volume per 1 gram of silica fume */
4997  heat_cf = 1000. / SYSIZE_POW3 * ( ( 1. / specgrav [ POZZ ] ) + ( float ) ( count [ POROSITY ] + count [ CH ] + count [ INERT ] ) / ( specgrav [ POZZ ] * ( float ) count [ POZZ ] ) ); //fixed
4998  }
4999 
5000  mass_fill_pozz = ( 1. - mass_agg ) * ( float ) ( count [ POZZ ] * specgrav [ POZZ ] ) / tot_mass;
5001  mass_fill = ( 1. - mass_agg ) * ( float ) ( count [ INERT ] * specgrav [ INERT ] +
5002  count [ ASG ] * specgrav [ ASG ] + count [ SLAG ] * specgrav [ SLAG ] +
5003  count [ CAS2 ] * specgrav [ CAS2 ] + count [ POZZ ] * specgrav [ POZZ ] +
5004  count [ CACL2 ] * specgrav [ CACL2 ] ) / tot_mass;
5005 #ifdef PRINTF
5006  printf("Calculated w/c is %.4f\n", w_to_c);
5007  printf("Calculated s/c is %.4f \n", s_to_c);
5008  printf("Calculated heat conversion factor is %f \n", heat_cf);
5009  printf("Calculated mass fractions of water and filler are %.4f and %.4f \n",
5011 #endif
5012  }
5013 
5014  molesdh2o = 0.0;
5015  alpha = 0.0; /* degree of hydration of clinker minerals*/
5016  /* heat4 contains measured heat release for C4AF hydration from */
5017  /* Fukuhara et al., Cem. and Conc. Res. article */
5018  heat4 = 0.0;
5019  mass_now = 0.0; /* total cement mass corrected for hydration */
5020  suminit = c3sinit + c2sinit + c3ainit + c4afinit;
5021  /* suminit=c3sinit+c2sinit+c3ainit+c4afinit+ncsbar+anhinit+heminit; */
5022  /* ctest is number of gypsum likely to form ettringite */
5023  /* 1 unit of C3A can react with 2.5 units of Gypsum */
5024  ctest = count [ DIFFGYP ];
5025 #ifdef PRINTF
5026  printf("ctest is %ld\n", ctest);
5027 #endif
5028  fflush(stdout);
5029  if ( ( float ) ctest > ( 2.5 * ( float ) ( count [ DIFFC3A ] + count [ DIFFC4A ] ) ) ) {
5030  ctest = ( long int ) ( 2.5 * ( float ) ( count [ DIFFC3A ] + count [ DIFFC4A ] ) );
5031  }
5032 
5033  for ( i = 0; i <= EMPTYP; i++ ) {
5034  if ( ( i != 0 ) && ( i <= ABSGYP ) && ( i != INERTAGG ) && ( i != CSH ) ) {
5035  heatsum += ( float ) count [ i ] * heatf [ i ] / molarv [ i ];
5036  /* Tabulate moles of H2O consumed by reactions so far */
5037  molesh2o += ( float ) count [ i ] * waterc [ i ] / molarv [ i ];
5038  }
5039 
5040  /* assume that all C3A which can, does form ettringite */
5041  if ( i == DIFFC3A ) {
5042  heatsum += ( ( float ) count [ DIFFC3A ] - ( float ) ctest / 2.5 ) * heatf [ C3A ] / molarv [ C3A ];
5043  }
5044 
5045  /* assume that all C4A which can, does form ettringite */
5046  if ( i == DIFFC4A ) {
5047  heatsum += ( ( float ) count [ DIFFC4A ] - ( float ) ctest / 2.5 ) * heatf [ C4AF ] / molarv [ C4AF ];
5048  }
5049 
5050  /* assume all gypsum which can, does form ettringite */
5051  /* rest will remain as gypsum */
5052  if ( i == DIFFGYP ) {
5053  heatsum += ( float ) ( count [ DIFFGYP ] - ctest ) * heatf [ GYPSUM ] / molarv [ GYPSUM ];
5054  /* 3.3 is the molar expansion from GYPSUM to ETTR */
5055  heatsum += ( float ) ctest * 3.30 * heatf [ ETTR ] / molarv [ ETTR ];
5056  molesdh2o += ( float ) ctest * 3.30 * waterc [ ETTR ] / molarv [ ETTR ];
5057  } else if ( i == DIFFCH ) {
5058  heatsum += ( float ) count [ DIFFCH ] * heatf [ CH ] / molarv [ CH ];
5059  molesdh2o += ( float ) count [ DIFFCH ] * waterc [ CH ] / molarv [ CH ];
5060  } else if ( i == DIFFFH3 ) {
5061  heatsum += ( float ) count [ DIFFFH3 ] * heatf [ FH3 ] / molarv [ FH3 ];
5062  molesdh2o += ( float ) count [ DIFFFH3 ] * waterc [ FH3 ] / molarv [ FH3 ];
5063  } else if ( i == DIFFCSH ) {
5064  /* use current CSH properties */
5065  heatsum += ( float ) count [ DIFFCSH ] * heatf [ CSH ] / molarvcsh [ cycle ];
5066  molesdh2o += ( float ) count [ DIFFCSH ] * watercsh [ cycle ] / molarvcsh [ cycle ];
5067  } else if ( i == DIFFETTR ) {
5068  heatsum += ( float ) count [ DIFFETTR ] * heatf [ ETTR ] / molarv [ ETTR ];
5069  molesdh2o += ( float ) count [ DIFFETTR ] * waterc [ ETTR ] / molarv [ ETTR ];
5070  } else if ( i == DIFFCACL2 ) {
5071  heatsum += ( float ) count [ DIFFCACL2 ] * heatf [ CACL2 ] / molarv [ CACL2 ];
5072  molesdh2o += ( float ) count [ DIFFCACL2 ] * waterc [ CACL2 ] / molarv [ CACL2 ];
5073  } else if ( i == DIFFAS ) {
5074  heatsum += ( float ) count [ DIFFAS ] * heatf [ ASG ] / molarv [ ASG ];
5075  molesdh2o += ( float ) count [ DIFFAS ] * waterc [ ASG ] / molarv [ ASG ];
5076  } else if ( i == DIFFCAS2 ) {
5077  heatsum += ( float ) count [ DIFFCAS2 ] * heatf [ CAS2 ] / molarv [ CAS2 ];
5078  molesdh2o += ( float ) count [ DIFFCAS2 ] * waterc [ CAS2 ] / molarv [ CAS2 ];
5079  }
5080  /* assume that all diffusing anhydrite leads to gypsum formation */
5081  else if ( i == DIFFANH ) {
5082  heatsum += ( float ) count [ DIFFANH ] * heatf [ GYPSUMS ] / molarv [ GYPSUMS ];
5083  /* 2 moles of water per mole of gypsum formed */
5084  molesdh2o += ( float ) count [ DIFFANH ] * 2.0 / molarv [ GYPSUMS ];
5085  }
5086  /* assume that all diffusing hemihydrate leads to gypsum formation */
5087  else if ( i == DIFFHEM ) {
5088  heatsum += ( float ) count [ DIFFHEM ] * heatf [ GYPSUMS ] / molarv [ GYPSUMS ];
5089  /* 1.5 moles of water per mole of gypsum formed */
5090  molesdh2o += ( float ) count [ DIFFHEM ] * 1.5 / molarv [ GYPSUMS ];
5091  } else if ( i == C3S ) {
5092  alpha += ( float ) ( c3sinit - count [ C3S ] );
5093  mass_now += specgrav [ C3S ] * ( float ) count [ C3S ];
5094  heat4 += .517 * ( float ) ( c3sinit - count [ C3S ] ) * specgrav [ C3S ];
5095  } else if ( i == C2S ) {
5096  alpha += ( float ) ( c2sinit - count [ C2S ] );
5097  mass_now += specgrav [ C2S ] * ( float ) count [ C2S ];
5098  heat4 += .262 * ( float ) ( c2sinit - count [ C2S ] ) * specgrav [ C2S ];
5099  } else if ( i == C3A ) {
5100  alpha += ( float ) ( c3ainit - count [ C3A ] );
5101  mass_now += specgrav [ C3A ] * ( float ) count [ C3A ];
5102  mc3ar = ( c3ainit - ( float ) count [ C3A ] ) / molarv [ C3A ];
5103  mc4ar = ( c4afinit - ( float ) count [ C4AF ] ) / molarv [ C4AF ];
5104  if ( ( mc3ar + mc4ar ) > 0.0 ) {
5105  frhyg = ( mc3ar / ( mc3ar + mc4ar ) ) * ( float ) count [ C3AH6 ] / molarv [ C3AH6 ];
5106  } else {
5107  frhyg = 0.0;
5108  }
5109 
5110  frettr = ( float ) count [ ETTR ] / molarv [ ETTR ];
5111  frafm = 3 * ( float ) count [ AFM ] / molarv [ AFM ];
5112  frtot = frafm + frettr + frhyg;
5113  if ( frtot > 0.0 ) {
5114  frettr /= frtot;
5115  frafm /= frtot;
5116  frhyg /= frtot;
5117  heat4 += frafm * 1.144 * ( float ) ( c3ainit - count [ C3A ] ) * specgrav [ C3A ];
5118  heat4 += frhyg * 0.908 * ( float ) ( c3ainit - count [ C3A ] ) * specgrav [ C3A ];
5119  heat4 += frettr * 1.672 * ( float ) ( c3ainit - count [ C3A ] ) * specgrav [ C3A ];
5120  }
5121  } else if ( i == C4AF ) {
5122  alpha += ( float ) ( c4afinit - count [ C4AF ] );
5123  mass_now += specgrav [ C4AF ] * ( float ) count [ C4AF ];
5124  mc3ar = ( c3ainit - ( float ) count [ C3A ] ) / molarv [ C3A ];
5125  mc4ar = ( c4afinit - ( float ) count [ C4AF ] ) / molarv [ C4AF ];
5126  if ( ( mc3ar + mc4ar ) > 0.0 ) {
5127  frhyg = ( mc4ar / ( mc3ar + mc4ar ) ) * ( float ) count [ C3AH6 ] / molarv [ C3AH6 ];
5128  } else {
5129  frhyg = 0.0;
5130  }
5131 
5132  frettr = ( float ) count [ ETTRC4AF ] / molarv [ ETTRC4AF ];
5133  frtot = frettr + frhyg;
5134  if ( frtot > 0.0 ) {
5135  frettr /= frtot;
5136  frhyg /= frtot;
5137  heat4 += frhyg * .418 * ( float ) ( c4afinit - count [ C4AF ] ) * specgrav [ C4AF ];
5138  heat4 += frettr * .725 * ( float ) ( c4afinit - count [ C4AF ] ) * specgrav [ C4AF ];
5139  }
5140  }
5141  /* else if(i==GYPSUM){
5142  * alpha+=(float)(ncsbar-count[GYPSUM]);
5143  * mass_now+=specgrav[GYPSUM]*(float)count[GYPSUM];
5144  * } */
5145  /* 0.187 kJ/g anhydrite for anhydrite --> gypsum conversion */
5146  else if ( i == ANHYDRITE ) {
5147  /* alpha+=(float)(anhinit-count[ANHYDRITE]);
5148  * mass_now+=specgrav[ANHYDRITE]*(float)count[ANHYDRITE]; */
5149  heat4 += .187 * ( float ) ( anhinit - count [ ANHYDRITE ] ) * specgrav [ ANHYDRITE ];
5150  /* 2 moles of water consumed per mole of anhydrite reacted */
5151  molesh2o += ( float ) ( anhinit - count [ ANHYDRITE ] ) * 2.0 / molarv [ ANHYDRITE ];
5152  }
5153  /* 0.132 kJ/g hemihydrate for hemihydrate-->gypsum conversion */
5154  else if ( i == HEMIHYD ) {
5155  /* alpha+=(float)(heminit-count[HEMIHYD]);
5156  * mass_now+=specgrav[HEMIHYD]*(float)count[HEMIHYD]; */
5157  heat4 += .132 * ( float ) ( heminit - count [ HEMIHYD ] ) * specgrav [ HEMIHYD ];
5158  /* 1.5 moles of water consumed per mole of anhydrite reacted */
5159  molesh2o += ( float ) ( heminit - count [ HEMIHYD ] ) * 1.5 / molarv [ HEMIHYD ];
5160  }
5161  }
5162 
5163  mass_fa_now = specgrav [ ASG ] * ( float ) count [ ASG ];
5164  mass_fa_now += specgrav [ CAS2 ] * ( float ) count [ CAS2 ];
5165  mass_fa_now += specgrav [ POZZ ] * ( float ) count [ POZZ ];
5166  if ( suminit != 0 ) {
5167  alpha = alpha / ( float ) suminit;
5168  } else {
5169  alpha = 0.0;
5170  }
5171 
5172  /* Current degree of hydration on a mass basis */
5173  if ( cemmass != 0.0 ) {
5174  alpha_cur = 1.0 - ( mass_now / cemmass );
5175  } else {
5176  alpha_cur = 0.0;
5177  }
5178 
5179  if ( flyashmass != 0.0 ) {
5180  alpha_fa_cur = 1.0 - ( mass_fa_now / flyashmass );
5181  } else {
5182  alpha_fa_cur = 0.0;
5183  }
5184 
5185  h2oinit = ( float ) porinit / molarv [ POROSITY ];
5186 
5187  /* Assume 780 J/g S for pozzolanic reaction */
5188  /* Each unit of silica fume consumes 1.35 units of CH, */
5189  /* so divide npr by 1.35 to get silca fume which has reacted */
5190  heat4 += 0.78 * ( ( float ) npr / 1.35 ) * specgrav [ POZZ ];
5191 
5192  /* Assume 800 J/g S for slag reaction */
5193  /* Seems reasonable with measurements of Biernacki and Richardson */
5194  heat4 += 0.8 * ( ( float ) nslagr ) * specgrav [ SLAG ];
5195 
5196  /* Assume 800 J/g AS for stratlingite formation (DeLarrard) */
5197  /* Each unit of AS consumes 1.3267 units of CH, */
5198  /* so divide nasr by 1.3267 to get ASG which has reacted */
5199  heat4 += 0.80 * ( ( float ) nasr / 1.3267 ) * specgrav [ ASG ];
5200 
5201  /* Should be additional code here for heat release due to CAS2 to */
5202  /* stratlingite conversion, but data unavailable at this time */
5203 
5204  /* Adjust heat sum for water left in system */
5205  water_left = ( long int ) ( ( h2oinit - molesh2o ) * molarv [ POROSITY ] + 0.5 );
5206  countkeep = count [ POROSITY ];
5207  heatsum += ( h2oinit - molesh2o - molesdh2o ) * heatf [ POROSITY ];
5208  if ( cyccnt == 0 ) {
5209 #ifdef OUTFILES
5210  fprintf(heatfile, "Cycle time(h) alpha_vol alpha_mass heat4(kJ/kg_solid) Gsratio2 G-s_ratio\n");
5211 #endif
5212  }
5213 
5214  heat_new = heat4; /* use heat4 for all adiabatic calculations */
5215  /* due to best agreement with calorimetry data */
5216 
5217  if ( cyccnt == 0 ) {
5218  //fprintf(chsfile,"Cycle time(h) alpha_mass Chemical shrinkage (ml/g cement)\n");
5219  }
5220 
5221  chs_new = ( ( float ) ( count [ EMPTYP ] + count [ POROSITY ] - water_left ) * heat_cf / 1000. );
5222  /* if((molesh2o>h2oinit)&&(sealed==1)){ */
5223  if ( ( ( water_left + water_off ) < 0 ) && ( sealed == 1 ) ) {
5224 #ifdef PRINTF
5225  printf("All water consumed at cycle %d \n", cyccnt);
5226  fflush(stdout);
5227 #endif
5228  //exit(1);the simulation can be stopped
5229  }
5230 
5231  //remove defined amount of water (outer drying)
5232  //makeinert(10);
5233 
5234  /* Attempt to create empty porosity to account for self-desiccation */
5235  if ( ( sealed == 1 ) && ( ( count [ POROSITY ] - water_left ) > 0 ) ) {
5236  poretodo = ( count [ POROSITY ] - pore_off ) - ( water_left - water_off );
5237  poretodo -= slagemptyp;
5238  if ( poretodo > 0 ) {
5240  poregone += poretodo;
5241  }
5242  }
5243 
5244  /* Output phase counts */
5245  /* phasfile for reactant and product phases */
5246 #ifdef OUTFILES
5247  if ( cyccnt == 0 ) {
5248  fprintf(phasfile, "#Cycle DoH Time Porosity C3S C2S C3A C4AF GYPSUM HEMIHYD ANHYDRITE POZZ INERT SLAG ASG CAS2 CH CSH C3AH6 ETTR ETTRC4AF AFM FH3 POZZCSH SLAGCSH CACL2 FREIDEL STRAT GYPSUMS CACO3 AFMC AGG ABSGYP EMPTYP HDCSH water_left \n");
5249  fprintf(perc_phases, "#Cyc Time[h] DoH SOL_per SOL_tot| Porosity C3S C2S C3A C4AF GYPSUM HEMIHYD ANHYDRITE POZZ INERT SLAG ASG CAS2 CH CSH C3AH6 ETTR ETTRC4AF AFM FH3 POZZCSH SLAGCSH CACL2 FREIDEL STRAT GYPSUMS CACO3 AFMC AGG ABSGYP EMPTYP HDCSH\n");
5250  }
5251 
5252  fprintf(phasfile, "%d %.4f %.4f ", cyccnt, alpha_cur, time_cur);
5253  fprintf(disprobfile, "%d %.4f %.4f ", cyccnt, alpha_cur, time_cur);
5254 
5255 
5257 
5258 
5259  for ( i = 0; i <= HDCSH; i++ ) {
5260  if ( ( i < DIFFCSH ) || ( i >= EMPTYP ) ) {
5261  if ( i == CSH ) {
5262  fprintf(phasfile, "%ld ", count [ CSH ] - count [ HDCSH ]);
5263  } else {
5264  fprintf(phasfile, "%ld ", count [ i ]);
5265  }
5266 
5267  fprintf(disprobfile, "%f ", disprob [ i ]);
5268  }
5269 
5270  // printf("%ld ",count[i]);
5271  }
5272 
5273  fprintf(disprobfile, "\n");
5274  fprintf(phasfile, "%ld\n", water_left);
5275  fflush(phasfile);
5276  fflush(disprobfile);
5277 #endif
5278 
5279 #ifdef PRINTF
5280  printf("\n");
5281 #endif
5282 
5283  if ( cycle == 0 ) {
5284  return;
5285  }
5286 
5287 #ifdef PRINTF
5288  #ifdef __TM_MODULE
5289  printf("Cycle cyccnt %d, DoH %f, time %f, GP %p\n", cyccnt, alpha_cur, time_cur, this->gp);
5290  #else
5291  printf("Cycle cyccnt %d, DoH %f, time %f\n", cyccnt, alpha_cur, time_cur);
5292  #endif
5293 #endif
5294 
5295  cyccnt += 1;
5296  fflush(stdout);
5297  /* Update current volume count for CH */
5298  chold = chnew;
5299  chnew = count [ CH ];
5300 
5301  /* See if ettringite is soluble yet */
5302  /* Gypsum 80% consumed, changed 06.09.00 from 90% to 80% */
5303  /* Gypsum 75% consumed, changed 09.09.01 from 80% to 75% */
5304  /* or system temperature exceeds 70 C */
5305  if ( ( ( ncsbar + anhinit + heminit ) != 0.0 ) || ( temp_cur >= 70.0 ) ) {
5306  /* Account for all sulfate sources and forms */
5307  if ( ( soluble [ ETTR ] == 0 ) && ( ( temp_cur >= 70.0 ) || ( count [ AFM ] != 0 ) ||
5308  ( ( ( ( float ) count [ GYPSUM ] + 1.42 * ( float ) count [ ANHYDRITE ] + 1.4 *
5309  ( float ) count [ HEMIHYD ] + ( float ) count [ GYPSUMS ] ) / ( ( float ) ncsbar +
5310  1.42 * ( float ) anhinit + 1.4 * ( float ) heminit + ( ( float ) netbar / 3.30 ) ) ) < 0.25 ) ) ) {
5311  soluble [ ETTR ] = 1;
5312 #ifdef PRINTF
5313  printf("Ettringite is soluble beginning at cycle %d \n", cycle);
5314 #endif
5315  /* identify all new soluble ettringite */
5316  passone(ETTR, ETTR, 2, 0);
5317  }
5318  }
5319 
5320  /* end of soluble ettringite test */
5321 
5322  /* Adjust ettringite solubility */
5323  /* if too many ettringites already in solution */
5324  if ( count [ DIFFETTR ] > DETTRMAX ) {
5325  disprob [ ETTR ] = 0.0;
5326  } else {
5327  disprob [ ETTR ] = disbase [ ETTR ];
5328  }
5329 
5330  /* Adjust CaCl2 solubility */
5331  /* if too many CaCl2 already in solution */
5332  if ( count [ DIFFCACL2 ] > DCACL2MAX ) {
5333  disprob [ CACL2 ] = 0.0;
5334  } else {
5335  disprob [ CACL2 ] = disbase [ CACL2 ];
5336  }
5337 
5338  /* Adjust CaCO3 solubility */
5339  /* if too many CaCO3 already in solution */
5340  if ( ( count [ DIFFCACO3 ] > DCACO3MAX ) && ( soluble [ ETTR ] == 0 ) ) {
5341  disprob [ CACO3 ] = 0.0;
5342  } else if ( count [ DIFFCACO3 ] > ( 4 * DCACO3MAX ) ) {
5343  disprob [ CACO3 ] = 0.0;
5344  } else {
5345  disprob [ CACO3 ] = disbase [ CACO3 ];
5346  }
5347 
5348  /* Adjust solubility of CH */
5349  /* based on amount of CH currently diffusing */
5350  /* Note that CH is always soluble to allow some */
5351  /* Ostwald ripening of the CH crystals */
5352  if ( ( float ) count [