Merge lp:~esys-p-dev/esys-particle/esys-darcy into lp:esys-particle

[Dion Weatherley]

Merge of esys-darcy into main trunk

Merging the integrated Darcy flow solver developed by Dr. Qi Shao with the main trunk of ESyS-Particle.

This merge incorporates Dr. Qi Shao's integrated Darcy flow solver with the main ESyS-Particle trunk

[Dion Weatherley (d-weatherley)]

review: Approve

=== added directory 'Doc/ESyS-Darcy'
=== added directory 'Doc/ESyS-Darcy/examples'
=== added directory 'Doc/ESyS-Darcy/examples/seepage'
=== added file 'Doc/ESyS-Darcy/examples/seepage/input.geo'
--- Doc/ESyS-Darcy/examples/seepage/input.geo	1970-01-01 00:00:00 +0000
+++ Doc/ESyS-Darcy/examples/seepage/input.geo	2019-07-28 14:55:57 +0000
@@ -0,0 +1,522 @@
+LSMGeometry 1.2
+BoundingBox 0 0 -0.01 0.04 0.02 0.07
+PeriodicBoundaries 0 0 0
+Dimension 3D
+BeginParticles
+Simple
+511
+0.00740403 0.00723698 0.00744442 0.00312945 0 0
+0.0171073 0.00808333 0.0162018 0.00331158 1 0
+0.00765378 0.00853487 0.0261773 0.00175397 2 0
+0.0159496 0.00777501 0.0349267 0.00332675 3 0
+0.00884835 0.00779284 0.0465758 0.00276472 4 0
+0.015543 0.00832276 0.0536022 0.00169508 5 0
+0.0221331 0.003481 0.0153119 0.00348576 6 0
+0.00596406 0.00553242 0.0407641 0.00150957 7 0
+0.0323902 0.00701856 0.0392024 0.00298129 8 0
+0.0150359 0.00834329 0.0453691 0.00234418 9 0
+0.00855937 0.00342381 0.0391418 0.00210383 10 0
+0.00179315 0.0182046 0.0142976 0.00179823 11 0
+0.0111781 0.00192011 0.0240379 0.00192042 12 0
+0.0196628 0.00689996 0.0112512 0.00241598 13 0
+0.024942 0.00842001 0.0262652 0.00308795 14 0
+0.027295 0.013957 0.016231 0.00150564 15 0
+0.00742343 0.0141679 0.0380271 0.00278792 16 0
+0.0216727 0.0102069 0.00749932 0.00304783 17 0
+0.0120856 0.0150207 0.0367895 0.00153555 18 0
+0.0100415 0.0168587 0.0426142 0.00314168 19 0
+0.00596843 0.0030228 0.0347412 0.00302346 20 0
+0.0112935 0.0031166 0.010052 0.00311947 21 0
+0.00751719 0.00466006 0.0432232 0.00152708 22 0
+0.0201785 0.0165272 0.0365523 0.00347492 23 0
+0.0214643 0.0106991 0.0341478 0.0029656 24 0
+0.0114554 0.0134664 0.0538252 0.00285964 25 0
+0.0314729 0.00639405 0.0446898 0.00262311 26 0
+0.0230855 0.00459689 0.022573 0.00255474 27 0
+0.0320717 0.0129698 0.0308292 0.00221269 28 0
+0.0230772 0.0168748 0.0194676 0.00274072 29 0
+0.0271645 0.00728577 0.013058 0.00255968 30 0
+0.0299018 0.0128437 0.0404557 0.00347676 31 0
+0.00332364 0.0151418 0.0276371 0.00332515 32 0
+0.0373931 0.00810864 0.0121803 0.00261654 33 0
+0.0332908 0.00940968 0.00967425 0.00242495 34 0
+0.0202158 0.0172934 0.00537452 0.00270846 35 0
+0.0371413 0.0066962 0.0458402 0.00285904 36 0
+0.0134183 0.00679717 0.00534724 0.00321175 37 0
+0.0382342 0.0164371 0.0259019 0.00176599 38 0
+0.0254657 0.0113559 0.0516564 0.00237733 39 0
+0.0200455 0.00177476 0.0198161 0.001776 40 0
+0.00302217 0.00803457 0.0245716 0.00302484 41 0
+0.0103181 0.00952543 0.0124983 0.00314521 42 0
+0.0384043 0.00722139 0.0515612 0.00159589 43 0
+0.0198441 0.00613669 0.0272463 0.00258811 44 0
+0.0349719 0.017209 0.0195774 0.00218768 45 0
+0.0147406 0.0180795 0.0440457 0.0019206 46 0
+0.0163462 0.0024843 0.0142653 0.00248488 47 0
+0.00448848 0.00247625 0.00782811 0.0024806 48 0
+0.00269461 0.0176299 0.0350801 0.00237049 49 0
+0.0162774 0.0115099 0.0127853 0.00160631 50 0
+0.0291583 0.00175865 0.0343068 0.00175969 51 0
+0.0176758 0.0117416 0.0384466 0.002257 52 0
+0.0113942 0.0144846 0.0332913 0.00207561 53 0
+0.0105094 0.00663566 0.0419762 0.00227142 54 0
+0.0277693 0.00849944 0.0207598 0.00313154 55 0
+0.00225796 0.0139071 0.0449231 0.00225884 56 0
+0.0327355 0.00869802 0.0524756 0.00210723 57 0
+0.0116614 0.0152128 0.0267079 0.00171612 58 0
+0.0128076 0.0117423 0.0356005 0.00178029 59 0
+0.0231826 0.00339772 0.0490703 0.00339785 60 0
+0.0257734 0.0132393 0.0362276 0.00245328 61 0
+0.0100857 0.018045 0.0282694 0.00188899 62 0
+0.00649234 0.0019779 0.0148566 0.001979 63 0
+0.0190784 0.0168062 0.0129325 0.00244603 64 0
+0.0274388 0.00283526 0.0186635 0.00283635 65 0
+0.00980859 0.011643 0.0238988 0.00266876 66 0
+0.00530997 0.0104109 0.0279439 0.00169793 67 0
+0.0150089 0.00230255 0.0188572 0.00230241 68 0
+0.0104202 0.00988296 0.034132 0.00158411 69 0
+0.0326759 0.0148408 0.036397 0.00183511 70 0
+0.00537875 0.017851 0.0127402 0.00214993 71 0
+0.0141031 0.0131392 0.0156461 0.00234038 72 0
+0.022045 0.0130303 0.0256106 0.00172639 73 0
+0.030883 0.0055074 0.0515429 0.00170148 74 0
+0.0171122 0.0156672 0.0189835 0.00282233 75 0
+0.0248191 0.0106904 0.0565557 0.0019531 76 0
+0.00998799 0.00230269 0.0303509 0.00230314 77 0
+0.0233557 0.00776256 0.0308694 0.0018403 78 0
+0.0153734 0.00613475 0.0407651 0.00277294 79 0
+0.0112266 0.0154402 0.00669456 0.0021556 80 0
+0.0284372 0.00330752 0.0246966 0.00330897 81 0
+0.0296624 0.0180387 0.0063879 0.00196153 82 0
+0.025597 0.0109553 0.00522174 0.0015032 83 0
+0.0115065 0.00724476 0.05159 0.00293994 84 0
+0.0266011 0.00629966 0.0541349 0.0033685 85 0
+0.0351079 0.00773157 0.0342485 0.00272183 86 0
+0.0229427 0.0116912 0.0474661 0.00252904 87 0
+0.00581224 0.00536645 0.0292007 0.0030036 88 0
+0.0305311 0.0104105 0.0564955 0.00278929 89 0
+0.0173384 0.00953158 0.0299498 0.00213646 90 0
+0.00391293 0.0139655 0.0328758 0.00207938 91 0
+0.00621244 0.0135989 0.0427517 0.00189021 92 0
+0.024043 0.00728268 0.0388119 0.00337211 93 0
+0.0186136 0.00722557 0.0482368 0.00237387 94 0
+0.0340442 0.00251141 0.0228404 0.00251184 95 0
+0.00491825 0.0183503 0.00613648 0.00165301 96 0
+0.0147637 0.012986 0.00900628 0.00273441 97 0
+0.0332292 0.00229225 0.0467351 0.00229396 98 0
+0.0282865 0.00300698 0.0408854 0.00300894 99 0
+0.00428008 0.00479348 0.0121972 0.00248558 100 0
+0.0083076 0.00619411 0.0218279 0.00322861 101 0
+0.0147392 0.00271629 0.0274105 0.00271994 102 0
+0.0376745 0.0177152 0.0565854 0.00228494 103 0
+0.019959 0.0183052 0.0166754 0.0016982 104 0
+0.038242 0.00860656 0.031193 0.00176244 105 0
+0.0242836 0.0151937 0.0399219 0.00198545 106 0
+0.0255509 0.0181291 0.0369613 0.00187095 107 0
+0.00838219 0.00505231 0.0135851 0.00185281 108 0
+0.0277088 0.00151638 0.049422 0.00151663 109 0
+0.0149918 0.00738756 0.0106314 0.00234305 110 0
+0.0365618 0.00566455 0.00834941 0.00200565 111 0
+0.01974 0.0115456 0.011923 0.00196263 112 0
+0.00685994 0.017412 0.0317913 0.00258912 113 0
+0.0104388 0.00590398 0.0334389 0.00245599 114 0
+0.0233839 0.00156916 0.0198111 0.0015725 115 0
+0.00785617 0.00833543 0.0376739 0.00307195 116 0
+0.0144027 0.0145482 0.0415487 0.00190752 117 0
+0.00974302 0.0178158 0.0355729 0.0021843 118 0
+0.0357762 0.0178467 0.052576 0.0021534 119 0
+0.00317369 0.00881829 0.0335775 0.00317579 120 0
+0.0135396 0.0150709 0.0229688 0.0024835 121 0
+0.0250317 0.0169316 0.032195 0.00307149 122 0
+0.0249193 0.0117534 0.0409116 0.00164436 123 0
+0.00230726 0.0131174 0.0368666 0.00230746 124 0
+0.0271016 0.0112457 0.0312579 0.00307011 125 0
+0.0338569 0.00788284 0.0245417 0.00294891 126 0
+0.015087 0.017037 0.0291434 0.00286593 127 0
+0.0235944 0.0133136 0.00457506 0.00165768 128 0
+0.021186 0.00823015 0.0440303 0.00266366 129 0
+0.0117129 0.00238318 0.00238331 0.00238526 130 0
+0.0106024 0.0112625 0.0075197 0.0020241 131 0
+0.0204841 0.00993927 0.0251963 0.00176161 132 0
+0.0269842 0.0122714 0.0107032 0.00244944 133 0
+0.0273927 0.00236758 0.0130094 0.00236832 134 0
+0.00319804 0.00965197 0.0416454 0.00319901 135 0
+0.0326107 0.00239774 0.0181546 0.0023979 136 0
+0.00187391 0.00652176 0.0378804 0.00187417 137 0
+0.0174243 0.00813536 0.0246287 0.001506 138 0
+0.0313724 0.00744182 0.0159657 0.00193022 139 0
+0.0263637 0.0133749 0.0196433 0.00206739 140 0
+0.0368989 0.00309771 0.0148131 0.00310353 141 0
+0.0331317 0.0184727 0.016621 0.00152978 142 0
+0.030803 0.0179462 0.0371897 0.00187845 143 0
+0.0384495 0.0155004 0.0356344 0.00155192 144 0
+0.030126 0.016506 0.0142283 0.00279956 145 0
+0.00464258 0.00791393 0.0149828 0.00172216 146 0
+0.0192921 0.00307373 0.0376961 0.00307455 147 0
+0.0305916 0.0068819 0.0292489 0.0028724 148 0
+0.0126891 0.0136203 0.0449236 0.00163748 149 0
+0.0105185 0.0177893 0.051347 0.00221122 150 0
+0.00318776 0.015905 0.0184879 0.00318971 151 0
+0.030311 0.0127497 0.0346018 0.00181544 152 0
+0.0296683 0.00533651 0.0336603 0.0018926 153 0
+0.00248031 0.0141469 0.0131478 0.00248124 154 0
+0.00216428 0.00952655 0.00619678 0.00216537 155 0
+0.0199893 0.00735672 0.0314117 0.00160125 156 0
+0.0268247 0.0183491 0.0212189 0.001651 157 0
+0.0222275 0.00288889 0.0428809 0.00288937 158 0
+0.00756801 0.0168824 0.0229114 0.00311793 159 0
+0.0312656 0.00623918 0.00333778 0.00333877 160 0
+0.0127276 0.0180344 0.0384488 0.00196729 161 0
+0.00185072 0.010301 0.028726 0.00185086 162 0
+0.0192515 0.0170171 0.0240123 0.00282793 163 0
+0.0088248 0.00553153 0.0170241 0.00164829 164 0
+0.0174593 0.0135792 0.0461189 0.00328253 165 0
+0.00808536 0.0140608 0.0190272 0.00171072 166 0
+0.0062331 0.0165946 0.0479246 0.00340578 167 0
+0.0357472 0.00829438 0.0419716 0.00155678 168 0
+0.0173856 0.0100909 0.0422682 0.00192304 169 0
+0.0230821 0.0140217 0.0102382 0.00185391 170 0
+0.00184968 0.0067189 0.0090493 0.00184971 171 0
+0.028432 0.0017601 0.0526059 0.00176021 172 0
+0.030613 0.00179218 0.0299505 0.00178407 173 0
+0.036528 0.0119222 0.0385787 0.00347314 174 0
+0.0185294 0.00300098 0.0316728 0.00300401 175 0
+0.022188 0.0175745 0.0429847 0.00242628 176 0
+0.020683 0.0130888 0.0515288 0.00232643 177 0
+0.0258479 0.00328883 0.0308183 0.00329078 178 0
+0.0234504 0.0090752 0.00245722 0.00247003 179 0
+0.00879877 0.0129209 0.0294136 0.00284995 180 0
+0.0258572 0.0142383 0.00195802 0.00196253 181 0
+0.0029827 0.00429277 0.0174766 0.00297787 182 0
+0.0256372 0.00738693 0.0338708 0.00182949 183 0
+0.0237353 0.00563493 0.00917512 0.00224192 184 0
+0.0251589 0.00158402 0.0532679 0.00158404 185 0
+0.00155557 0.00526344 0.042439 0.00155593 186 0
+0.0354467 0.00743304 0.0501311 0.00169709 187 0
+0.022148 0.0170426 0.0483429 0.00295866 188 0
+0.00161019 0.00161185 0.0105829 0.00161257 189 0
+0.0166357 0.0034659 0.0522416 0.00346604 190 0
+0.0280514 0.0176819 0.0278071 0.00231822 191 0
+0.038291 0.00872478 0.00155462 0.001551 192 0
+0.00652405 0.0033354 0.00291777 0.00292349 193 0
+0.0148099 0.00168521 0.0387253 0.00168542 194 0
+0.0325501 0.0177729 0.0274042 0.00222942 195 0
+0.00234334 0.00234271 0.0486866 0.0023434 196 0
+0.0235571 0.0178086 0.00889742 0.00219226 197 0
+0.00650385 0.00776521 0.0176745 0.00158606 198 0
+0.038139 0.018139 0.014236 0.00186139 199 0
+0.00560776 0.00640613 0.0529907 0.00200518 200 0
+0.0327698 0.00168103 0.0412371 0.00168104 201 0
+0.00321623 0.00782948 0.0485721 0.00321733 202 0
+0.00168006 0.00167935 0.0138615 0.00168074 203 0
+0.0120099 0.0122088 0.00289867 0.00290008 204 0
+0.00684948 0.0184589 0.0272349 0.0015438 205 0
+0.0172233 0.0124645 0.0560698 0.00341054 206 0
+0.0307111 0.0174269 0.0319002 0.00257369 207 0
+0.0317224 0.0131959 0.0056655 0.0033502 208 0
+0.0223864 0.0121828 0.0385944 0.00181601 209 0
+0.0333741 0.00162267 0.0379904 0.00162341 210 0
+0.0126385 0.00474775 0.0451379 0.00198505 211 0
+0.0333118 0.00975449 0.047728 0.00226902 212 0
+0.00175565 0.0182005 0.0390556 0.00175591 213 0
+0.00205201 0.00539685 0.00538364 0.00205289 214 0
+0.00323014 0.0146552 0.0406688 0.00189881 215 0
+0.0264494 0.0149921 0.0068019 0.00233734 216 0
+0.0273275 0.0180763 0.0176918 0.00192378 217 0
+0.0240693 0.00191583 0.0385519 0.00191617 218 0
+0.00916667 0.0124688 0.0453142 0.0020909 219 0
+0.0118018 0.00792646 0.0179651 0.00228676 220 0
+0.00211426 0.0020685 0.0302846 0.00206859 221 0
+0.0325379 0.00517236 0.00846175 0.00205096 222 0
+0.0174031 0.00675963 0.0570005 0.00237341 223 0
+0.0308439 0.00335187 0.036816 0.00166135 224 0
+0.0378043 0.0117057 0.00526434 0.00219581 225 0
+0.0187597 0.0183295 0.0450964 0.00167053 226 0
+0.0306812 0.0113159 0.0255473 0.00183528 227 0
+0.0257224 0.014471 0.0443686 0.00247729 228 0
+0.00564528 0.0136085 0.00794232 0.00348379 229 0
+0.00690026 0.0016242 0.0110379 0.00162423 230 0
+0.010398 0.002126 0.0161074 0.0021265 231 0
+0.0291003 0.00998606 0.0491639 0.00218657 232 0
+0.01055 0.00175199 0.0343268 0.0017523 233 0
+0.0146657 0.0174941 0.0344697 0.00250615 234 0
+0.0367023 0.0129797 0.0230008 0.00300049 235 0
+0.0170752 0.00958292 0.00687754 0.00163407 236 0
+0.0305487 0.0140316 0.0103414 0.00154308 237 0
+0.012143 0.00512642 0.037407 0.00193262 238 0
+0.0144416 0.00239774 0.0348906 0.00223713 239 0
+0.0208874 0.0152295 0.0161386 0.00155911 240 0
+0.017102 0.0119374 0.023918 0.00239269 241 0
+0.0380168 0.0116568 0.0333423 0.0019839 242 0
+0.020783 0.017238 0.0536266 0.00232437 243 0
+0.00459883 0.00936175 0.0108443 0.00177585 244 0
+0.00979842 0.011681 0.04944 0.0021586 245 0
+0.0157955 0.0168473 0.0513997 0.00315337 246 0
+0.0222891 0.0142212 0.0562616 0.00195534 247 0
+0.0328359 0.00573027 0.0487517 0.00171545 248 0
+0.0233948 0.00708164 0.0189688 0.00181446 249 0
+0.0167859 0.00612698 0.0211226 0.00201806 250 0
+0.0366613 0.0134709 0.0108203 0.00298431 251 0
+0.0117147 0.0182589 0.0250377 0.00174352 252 0
+0.0351538 0.0178268 0.0238751 0.00217522 253 0
+0.0343496 0.0125793 0.0176437 0.00286872 254 0
+0.0311049 0.0052083 0.0205014 0.00156194 255 0
+0.0382786 0.0182785 0.0178056 0.00172253 256 0
+0.0280292 0.003005 0.0077153 0.00301014 257 0
+0.00156874 0.0110855 0.00953688 0.0015688 258 0
+0.0336558 0.0173626 0.040316 0.00240105 259 0
+0.0101468 0.0129838 0.0162281 0.00166167 260 0
+0.0365685 0.0163658 0.0311024 0.00343461 261 0
+0.0084464 0.018151 0.00767331 0.00184901 262 0
+0.00183102 0.00828557 0.0128648 0.00183638 263 0
+0.0144202 0.0171772 0.0580798 0.00242996 264 0
+0.037643 0.0101872 0.0274533 0.00235827 265 0
+0.0306139 0.00214476 0.0558196 0.00214483 266 0
+0.0369547 0.0125392 0.0450446 0.00304545 267 0
+0.0137442 0.0118058 0.0271828 0.00231246 268 0
+0.00572437 0.0112326 0.0138614 0.00195718 269 0
+0.013553 0.0174369 0.00257317 0.00256338 270 0
+0.0178381 0.0167656 0.0416959 0.00218456 271 0
+0.0327479 0.0174337 0.00964386 0.00256805 272 0
+0.0115658 0.0170655 0.018374 0.00293691 273 0
+0.00722468 0.00201345 0.0187804 0.00201356 274 0
+0.0380832 0.00191675 0.0263024 0.00191783 275 0
+0.00151268 0.0184873 0.0220424 0.00151274 276 0
+0.0188008 0.0134214 0.0281418 0.00240862 277 0
+0.0370407 0.0108358 0.053941 0.00278181 278 0
+0.0382324 0.0149963 0.0189764 0.00176843 279 0
+0.0333946 0.00268916 0.0333306 0.00269044 280 0
+0.0379133 0.00495567 0.0316106 0.00192929 281 0
+0.0139653 0.00453043 0.031477 0.00181791 282 0
+0.0165215 0.00176899 0.0416556 0.00171339 283 0
+0.0233985 0.00402321 0.0266658 0.00159446 284 0
+0.0289804 0.0135746 0.047063 0.00184511 285 0
+0.0276874 0.00829591 0.0427048 0.00206229 286 0
+0.0104164 0.0050226 0.0267236 0.00225462 287 0
+0.0299703 0.0043149 0.0152598 0.00156971 288 0
+0.0355399 0.00830889 0.00553889 0.00189664 289 0
+0.0211319 0.00999705 0.0288925 0.0018021 290 0
+0.012526 0.0183266 0.0546273 0.00167363 291 0
+0.00215674 0.0178438 0.0429304 0.00215679 292 0
+0.036927 0.0169288 0.00589332 0.00307357 293 0
+0.00711901 0.0106929 0.0199971 0.00177195 294 0
+0.00608676 0.0124662 0.0516043 0.00212781 295 0
+0.0332294 0.00995106 0.0430838 0.0016568 296 0
+0.00160766 0.00453795 0.0350303 0.00160967 297 0
+0.0326308 0.0131346 0.0510869 0.00254587 298 0
+0.0234443 0.0110963 0.0157333 0.00331821 299 0
+0.0161996 0.00450485 0.0462712 0.00176098 300 0
+0.0217518 0.00172445 0.0286186 0.00161838 301 0
+0.0309133 0.0112472 0.0136353 0.00255638 302 0
+0.00872013 0.0113484 0.0413458 0.00175616 303 0
+0.0352245 0.0110546 0.0303423 0.0015085 304 0
+0.037468 0.00814346 0.0204297 0.00253284 305 0
+0.0383311 0.0150886 0.0537058 0.00166904 306 0
+0.0120816 0.00233167 0.0416034 0.00233396 307 0
+0.0294322 0.00800664 0.025063 0.00150829 308 0
+0.0343668 0.0147665 0.0556005 0.00225467 309 0
+0.0207971 0.0176529 0.0289065 0.00235005 310 0
+0.0319273 0.00217282 0.0128965 0.00217288 311 0
+0.0371121 0.00416696 0.0037222 0.002892 312 0
+0.00775536 0.00188664 0.0263674 0.00188868 313 0
+0.0331224 0.0174383 0.0483866 0.00256238 314 0
+0.0112521 0.0181341 0.0318277 0.00186785 315 0
+0.0282315 0.0149082 0.0508115 0.00220593 316 0
+0.0228505 0.00329858 0.00329885 0.00329944 317 0
+0.0178559 0.00243378 0.0233342 0.00243497 318 0
+0.00164242 0.0183563 0.00638179 0.00164843 319 0
+0.0186533 0.00744749 0.00287275 0.00262119 320 0
+0.0347361 0.0180106 0.0359326 0.00199199 321 0
+0.0196083 0.00662829 0.0231802 0.00152232 322 0
+0.0334753 0.00183238 0.00183527 0.00183807 323 0
+0.0105506 0.0160335 0.0122247 0.00338778 324 0
+0.0287291 0.00917687 0.0360068 0.00234093 325 0
+0.0231773 0.0134731 0.0290722 0.001944 326 0
+0.0204786 0.0099285 0.0207758 0.00265931 327 0
+0.0259044 0.00788416 0.0467202 0.00235282 328 0
+0.0221521 0.00567684 0.0339741 0.00207081 329 0
+0.00220569 0.0130599 0.0495825 0.00220624 330 0
+0.0166981 0.00261089 0.00260843 0.0026122 331 0
+0.0352549 0.00744536 0.0295519 0.00184709 332 0
+0.0181356 0.00514363 0.00685589 0.00201028 333 0
+0.0279339 0.0105814 0.00227513 0.00227643 334 0
+0.0369969 0.0030031 0.0503873 0.0030034 335 0
+0.0206375 0.00897867 0.0569603 0.00154939 336 0
+0.0382062 0.00179401 0.00760958 0.0017943 337 0
+0.0271071 0.0149506 0.0558297 0.00293743 338 0
+0.0354024 0.00407603 0.0263096 0.00153291 339 0
+0.0187072 0.00154247 0.0481111 0.00154347 340 0
+0.0150173 0.0163305 0.00624998 0.00155225 341 0
+0.00686868 0.0102587 0.00260768 0.0026095 342 0
+0.0100899 0.0169072 0.0566386 0.00179012 343 0
+0.0137025 0.00746051 0.0286139 0.00226749 344 0
+0.0129312 0.0109103 0.0315169 0.00220639 345 0
+0.0158624 0.0142056 0.00150858 0.00150852 346 0
+0.017519 0.0183293 0.0556543 0.00167074 347 0
+0.0079414 0.00892021 0.031694 0.00183696 348 0
+0.0382482 0.0182475 0.0214964 0.00175278 349 0
+0.0259955 0.00656889 0.0169037 0.0015487 350 0
+0.0161249 0.00187418 0.00996327 0.00187507 351 0
+0.00865046 0.0101062 0.0169744 0.00166719 352 0
+0.006204 0.0121718 0.0171857 0.00153468 353 0
+0.0384895 0.0184891 0.0502027 0.00151113 354 0
+0.0161742 0.0184681 0.0388438 0.00153249 355 0
+0.0117072 0.00161091 0.0374811 0.00161099 356 0
+0.0197249 0.0126478 0.00269504 0.00270316 357 0
+0.0379822 0.010017 0.0493093 0.00201857 358 0
+0.0079389 0.00970356 0.0546289 0.00235297 359 0
+0.00634777 0.0147387 0.0148706 0.00174477 360 0
+0.0160297 0.0113212 0.0199518 0.00177307 361 0
+0.0262327 0.00692446 0.00516674 0.00200391 362 0
+0.0190803 0.00167763 0.0058822 0.00155 363 0
+0.021295 0.00167028 0.0534442 0.00167044 364 0
+0.0354161 0.0165134 0.0159228 0.00156673 365 0
+0.0259647 0.0183665 0.0505247 0.00163387 366 0
+0.00323668 0.00262533 0.0396654 0.00262566 367 0
+0.00170074 0.0182995 0.0311909 0.00170109 368 0
+0.0202307 0.00595664 0.0192722 0.00155746 369 0
+0.0380636 0.00591639 0.0275328 0.00193819 370 0
+0.0353452 0.00778029 0.0163646 0.00207745 371 0
+0.0223714 0.00840571 0.0502387 0.00180997 372 0
+0.00221456 0.00647535 0.054996 0.00193698 373 0
+0.0289452 0.00529415 0.0483547 0.00196326 374 0
+0.0287647 0.00832274 0.00847278 0.00242411 375 0
+0.0172068 0.0179303 0.00898301 0.0020752 376 0
+0.00202381 0.00924055 0.00202347 0.00202428 377 0
+0.0245742 0.0170822 0.0245145 0.00253671 378 0
+0.0379998 0.00200056 0.0392715 0.00200151 379 0
+0.0334767 0.00593786 0.0121811 0.00186137 380 0
+0.0180743 0.00181884 0.0448255 0.00182027 381 0
+0.0122542 0.0108769 0.0466382 0.00162706 382 0
+0.00655569 0.0168082 0.0539671 0.00264169 383 0
+0.0383407 0.0183407 0.00165259 0.00166263 384 0
+0.0214504 0.00172539 0.0252907 0.00172862 385 0
+0.00214564 0.0175114 0.00999983 0.00215176 386 0
+0.0150882 0.0183661 0.0248541 0.00163821 387 0
+0.0384032 0.00870137 0.00750029 0.00159694 388 0
+0.0161057 0.010794 0.00405664 0.00158625 389 0
+0.0328496 0.00829805 0.0196323 0.00209146 390 0
+0.0156947 0.00564037 0.0244717 0.00154824 391 0
+0.00195455 0.00195538 0.00195608 0.00195661 392 0
+0.0372154 0.0145755 0.0499549 0.00227933 393 0
+0.0243971 0.0170321 0.0139171 0.002969 394 0
+0.00762674 0.00305111 0.0497952 0.00304934 395 0
+0.0379301 0.00482145 0.0235776 0.00207029 396 0
+0.0264092 0.00350722 0.035982 0.00190991 397 0
+0.0193412 0.0144205 0.00912034 0.00207234 398 0
+0.035643 0.00411888 0.0370425 0.00188387 399 0
+0.0297379 0.0104753 0.0452151 0.00184417 400 0
+0.00583399 0.00269151 0.0453085 0.00180035 401 0
+0.0358471 0.00253241 0.0293966 0.00195015 402 0
+0.0182762 0.0131768 0.0149761 0.00179601 403 0
+0.0350586 0.00545458 0.0562953 0.00341809 404 0
+0.00218552 0.0128416 0.0227735 0.00218788 405 0
+0.0152214 0.0178986 0.0150455 0.00209452 406 0
+0.0353873 0.0050777 0.0187091 0.00150126 407 0
+0.0384627 0.0015372 0.0314742 0.00153768 408 0
+0.00211531 0.00211565 0.044241 0.00211629 409 0
+0.0313858 0.0151228 0.0236641 0.00247037 410 0
+0.00923771 0.00179787 0.0445522 0.00179785 411 0
+0.0383583 0.0108847 0.0152595 0.00164174 412 0
+0.0238522 0.0178169 0.00218615 0.00218987 413 0
+0.014146 0.00867598 0.0210404 0.00165307 414 0
+0.037655 0.00201786 0.0347177 0.0018402 415 0
+0.0320816 0.00203275 0.050898 0.0020332 416 0
+0.0154327 0.0151909 0.0126162 0.00155675 417 0
+0.0313455 0.00965986 0.032855 0.00174282 418 0
+0.0281756 0.0177571 0.0475384 0.00213551 419 0
+0.0291461 0.0181753 0.0404917 0.00182473 420 0
+0.0186263 0.00174812 0.0567437 0.0017482 421 0
+0.0329754 0.0176709 0.00230144 0.00230164 422 0
+0.0127313 0.00532341 0.0145625 0.00212689 423 0
+0.0376574 0.0176587 0.0465346 0.00234295 424 0
+0.0271144 0.0128908 0.0239936 0.00236246 425 0
+0.0215385 0.00459648 0.0562059 0.00236118 426 0
+0.00738788 0.00904432 0.0507336 0.00163676 427 0
+0.0382684 0.00173155 0.0547663 0.00173191 428 0
+0.0257193 0.0183717 0.0408645 0.00162862 429 0
+0.0126253 0.00219756 0.04845 0.00219787 430 0
+0.0109563 0.00196688 0.0201579 0.00196749 431 0
+0.0254353 0.0184314 0.045143 0.00156921 432 0
+0.00570082 0.0109462 0.045671 0.00171668 433 0
+0.0282435 0.0110536 0.0168082 0.00160284 434 0
+0.00179531 0.00403575 0.0269695 0.00180045 435 0
+0.016429 0.0137832 0.00513348 0.00157774 436 0
+0.0126431 0.00782114 0.02426 0.00200227 437 0
+0.00241284 0.0105334 0.0534433 0.00241296 438 0
+0.00536789 0.0179582 0.0402004 0.00204203 439 0
+0.00729154 0.018056 0.0163168 0.00194866 440 0
+0.00413258 0.00214982 0.0535375 0.00214994 441 0
+0.0236329 0.00873546 0.0114316 0.00160723 442 0
+0.00264859 0.00264829 0.0228373 0.00264955 443 0
+0.0378351 0.00650969 0.0394615 0.00216522 444 0
+0.0380916 0.018083 0.0105314 0.0018534 445 0
+0.0321456 0.0132646 0.0452373 0.00182279 446 0
+0.00966394 0.00166826 0.00580348 0.00166644 447 0
+0.00168093 0.0046194 0.0519373 0.00168113 448 0
+0.0347692 0.0102218 0.00210724 0.00210747 449 0
+0.00255495 0.0172316 0.00251778 0.00252309 450 0
+0.0112111 0.00226283 0.0536347 0.0022629 451 0
+0.0129764 0.00955392 0.0551326 0.00153717 452 0
+0.0132403 0.00456442 0.0220498 0.00197873 453 0
+0.0299834 0.0143163 0.0183655 0.00188999 454 0
+0.032864 0.00190089 0.0270448 0.00190282 455 0
+0.0068479 0.00166642 0.0229382 0.00166951 456 0
+0.0203736 0.00242782 0.00956418 0.00243119 457 0
+0.00675756 0.00172359 0.0420532 0.00172428 458 0
+0.0311484 0.0119802 0.0208795 0.00173721 459 0
+0.0222478 0.00944784 0.053868 0.00196903 460 0
+0.0168253 0.01213 0.0508497 0.00170684 461 0
+0.0183685 0.00504764 0.0437566 0.00160412 462 0
+0.0232737 0.0131858 0.02223 0.00187716 463 0
+0.0293226 0.0140501 0.00154323 0.00155941 464 0
+0.0187692 0.0092528 0.0520583 0.0019576 465 0
+0.0365652 0.0138516 0.0019462 0.00194697 466 0
+0.0215012 0.0130162 0.0424346 0.00221918 467 0
+0.0381017 0.0144179 0.0153672 0.00190392 468 0
+0.00332102 0.0129332 0.00189674 0.00189745 469 0
+0.0177389 0.0182206 0.00177586 0.00178121 470 0
+0.0307391 0.017877 0.0200459 0.00212367 471 0
+0.0363481 0.00156879 0.0191933 0.00156882 472 0
+0.0290274 0.0184173 0.0235843 0.00158279 473 0
+0.0339247 0.0137638 0.0271061 0.00201849 474 0
+0.0283487 0.0176869 0.00233598 0.00231694 475 0
+0.0282852 0.0184908 0.0350055 0.00151042 476 0
+0.0330134 0.00185311 0.00549684 0.00185315 477 0
+0.0242935 0.00160266 0.010027 0.00160277 478 0
+0.030886 0.0172073 0.044091 0.00228919 479 0
+0.0312084 0.0175061 0.0533512 0.00249442 480 0
+0.0341718 0.0183149 0.0135427 0.00168831 481 0
+0.0120748 0.0105452 0.0429707 0.00205716 482 0
+0.012124 0.0120933 0.0196146 0.00222838 483 0
+0.0147474 0.0184644 0.0115119 0.00154335 484 0
+0.0378522 0.0178499 0.0386472 0.00215028 485 0
+0.00192368 0.0146758 0.0546279 0.00192372 486 0
+0.0155852 0.00901374 0.0502825 0.00169673 487 0
+0.00803869 0.0170089 0.00298788 0.0029914 488 0
+0.0362089 0.00233206 0.0431987 0.00233296 489 0
+0.0330852 0.010437 0.0281984 0.00158274 490 0
+0.0122583 0.0182774 0.0467047 0.00172264 491 0
+0.00226078 0.0177388 0.0517969 0.00226149 492 0
+0.00280152 0.00993607 0.018667 0.00280199 493 0
+0.0146218 0.00189707 0.00652101 0.00189122 494 0
+0.0384981 0.0184995 0.0351118 0.00150211 495 0
+0.00425058 0.00164772 0.0266976 0.00164977 496 0
+0.0133873 0.0182043 0.00848053 0.00179734 497 0
+0.0383755 0.00427334 0.0191478 0.00162463 498 0
+0.0383818 0.00897087 0.0576935 0.00161829 499 0
+0.0382429 0.0137606 0.0571964 0.00175716 500 0
+0.0168499 0.0133127 0.0326882 0.00254127 501 0
+0.0150799 0.00917761 0.00156015 0.00156038 502 0
+0.0246948 0.0182753 0.0536291 0.00172506 503 0
+0.0350383 0.0021399 0.00991198 0.00214006 504 0
+0.00725777 0.0119493 0.0338505 0.00194902 505 0
+0.0277642 0.0178239 0.0100572 0.00217661 506 0
+0.00168331 0.00557014 0.00168047 0.00168212 507 0
+0.011902 0.0124199 0.0392896 0.0020792 508 0
+0.00526411 0.012825 0.0557359 0.00190477 509 0
+0.0275611 0.00235525 0.0459032 0.00210683 510 0
+EndParticles
+BeginConnect
+0
+EndConnect

=== added file 'Doc/ESyS-Darcy/examples/seepage/run.sh'
--- Doc/ESyS-Darcy/examples/seepage/run.sh	1970-01-01 00:00:00 +0000
+++ Doc/ESyS-Darcy/examples/seepage/run.sh	2019-07-28 14:55:57 +0000
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+#a bash file to run the simulation and create vtk files for visualization in Paraview
+#use the command "sh run.sh" in terminal for execution 
+
+
+if [ -d "results" ]; then
+  rm -r results
+fi
+
+mkdir results
+
+mpirun -np 2 `which esysparticle` upward_seepage.py
+
+dump2vtk -i results/snap -o results/particles_ -rot -t 0 6 100000
+
+printf "\nSimulation done, results saved in /results\n"

=== added file 'Doc/ESyS-Darcy/examples/seepage/upward_seepage.py'
--- Doc/ESyS-Darcy/examples/seepage/upward_seepage.py	1970-01-01 00:00:00 +0000
+++ Doc/ESyS-Darcy/examples/seepage/upward_seepage.py	2019-07-28 14:55:57 +0000
@@ -0,0 +1,204 @@
+#upward_seepage.py: A upward seepage simulation through particle assembly using ESyS-Darcy (a fluid-coupled version of ESyS-Particle)
+#       Author: Q.Shao
+#       Date:  21 August 2017
+#       Organisation: University of Queensland
+#       (C) All rights reserved, 2017.
+#
+#
+				
+#import the appropriate ESyS-Particle modules:
+from esys.lsm import *
+from esys.lsm.util import *
+from esys.lsm.geometry import *
+
+#instantiate a simulation object
+#and initialise the neighbour search algorithm:
+sim = LsmMpi(numWorkerProcesses=1, mpiDimList=[1,1,1])
+
+sim.initNeighbourSearch(
+    particleType="RotSphere",
+    gridSpacing=0.008,
+    verletDist=0.0003
+)
+
+#specify the number of timesteps and timestep increment:
+sim.setNumTimeSteps(500000)
+sim.setTimeStepSize(1.0e-6)
+
+#specify the spatial domain for the simulation:
+domain = BoundingBox(Vec3(0.,0.,-0.01), Vec3(0.04,0.02,0.07))
+sim.setSpatialDomain(domain)
+
+#add fluid to the domain:
+sim.createFluidInteraction(
+  cellside=0.01,         #side length of fluid cell
+  Bw=1.0e6,              #bulk modulus of fluid
+  Bp=1.0e6,              #bulk modulus of particle
+  Mu=1.0e-3,             #fluid viscosity
+  alpha=0.5,             #adjusting factor between two time steps
+  flowrate=0.1,          #rate of inflow
+  pressure=0.,           #hydraulic pressure
+  inflow=Vec3(0,0,1),    #directions of inflows
+  outflow=Vec3(0,0,1),   #directions of outflows 
+  fluid_timestep=1.0e-4  #time step size for updating fluid phase
+)
+
+#read in geometry input file
+sim.readGeometry("input.geo")
+
+#set partilce density
+sim.setParticleDensity (
+  tag = 0,
+  mask = -1,
+  Density = 2600
+)
+
+#add a left side wall to the model:
+sim.createWall (
+  name = "left_wall",
+  posn = Vec3(0.0000, 0.0000, 0.0000),
+  normal = Vec3(1.0000, 0.0000, 0.0000)
+)
+
+#add a right side wall to the model:
+sim.createWall (
+  name = "right_wall",
+  posn = Vec3(0.04, 0.0000, 0.0000),
+  normal = Vec3(-1.0000, 0.0000, 0.0000)
+)
+
+#add a back side wall to the model:
+sim.createWall (
+  name = "back_wall",
+  posn = Vec3(0.0000, 0.0000, 0.0000),
+  normal = Vec3(0.0000, 1.0000, 0.0000)
+)
+
+#add a front side wall to the model:
+sim.createWall (
+  name = "front_wall",
+  posn = Vec3(0.0000, 0.02, 0.0000),
+  normal = Vec3(0.0000, -1.0000, 0.0000)
+)
+
+#add a bottom side wall to the model:
+sim.createWall (
+  name = "bottom_wall",
+  posn = Vec3(0.0000, 0.0000, 0.0000),
+  normal = Vec3(0.0000, 0.0000, 1.0000)
+)
+
+#specify that particles undergo frictional interactions:
+sim.createInteractionGroup (
+  RotFrictionPrms (
+    name = "friction",
+    normalK = 1.0e6,
+    dynamicMu = 0.4,
+    staticMu = 0.6,
+    shearK = 1.0e6,
+    scaling = True,
+    rigid = True,
+    meanR_scaling = True
+  )
+)
+
+#specify that particles undergo elastic repulsion
+#from the left side wall:
+sim.createInteractionGroup (
+  NRotElasticWallPrms (
+    name = "lw_repel",
+    wallName = "left_wall",
+    normalK = 1.0e6
+  )
+)
+
+#specify that particles undergo elastic repulsion
+#from the right side wall:
+sim.createInteractionGroup (
+  NRotElasticWallPrms (
+    name = "rw_repel",
+    wallName = "right_wall",
+    normalK = 1.0e6
+  )
+)
+
+#specify that particles undergo elastic repulsion
+#from the back side wall:
+sim.createInteractionGroup (
+  NRotElasticWallPrms (
+    name = "bkw_repel",
+    wallName = "back_wall",
+    normalK = 1.0e6
+  )
+)
+
+
+#specify that particles undergo elastic repulsion
+#from the front side wall:
+sim.createInteractionGroup (
+  NRotElasticWallPrms (
+    name = "fw_repel",
+    wallName = "front_wall",
+    normalK = 1.0e6
+  )
+)
+
+
+#specify that particles undergo elastic repulsion
+#from the bottom side wall:
+sim.createInteractionGroup (
+  NRotElasticWallPrms (
+    name = "btw_repel",
+    wallName = "bottom_wall",
+    normalK = 1.0e6
+  )
+)
+
+#add a CheckPointer to store simulation data
+sim.createCheckPointer (
+  CheckPointPrms (
+  fileNamePrefix = "results/snap",
+    beginTimeStep=0,
+    endTimeStep=500000,
+    timeStepIncr=100000
+  )
+)
+
+#create a FieldSaver to store distributed pore pressures
+sim.createFieldSaver (
+  FluidScalarFieldSaverPrms(
+    fieldName="disP",
+    fileName="results/Pressure",
+    fileFormat="VTI",
+    beginTimeStep=0,
+    endTimeStep=500000,
+    timeStepIncr=100000
+  )
+)
+
+#create a FieldSaver to store porosities
+sim.createFieldSaver (
+  FluidScalarFieldSaverPrms(
+    fieldName="Phi",
+    fileName="results/Porosity",
+    fileFormat="VTI",
+    beginTimeStep=0,
+    endTimeStep=500000,
+    timeStepIncr=100000
+  )
+)
+
+#create a FieldSaver to store permeabilities
+sim.createFieldSaver (
+  FluidScalarFieldSaverPrms(
+    fieldName="K",
+    fileName="results/Permeability",
+    fileFormat="VTI",
+    beginTimeStep=0,
+    endTimeStep=500000,
+    timeStepIncr=100000
+  )
+)
+
+#Execute the simulation:
+sim.run()

=== modified file 'Fields/FieldMaster.cpp'
--- Fields/FieldMaster.cpp	2017-01-04 08:26:43 +0000
+++ Fields/FieldMaster.cpp	2019-07-28 14:55:57 +0000
@@ -68,6 +68,18 @@
     m_write_type=WRITE_TYPE_RAW_WITH_ID;
   }else if(savetype=="RAW_WITH_POS_ID"){
     m_write_type=WRITE_TYPE_RAW_WITH_POS_ID;
+  /****fluid contents: begin****/
+  }else if(savetype=="RAW_WITH_POS"){
+    m_write_type=WRITE_TYPE_RAW_WITH_POS;
+  }else if(savetype=="RAW_WITH_ID_POS"){
+    m_write_type=WRITE_TYPE_RAW_WITH_ID_POS;
+  }else if(savetype=="RAW_WITH_PARTICLE"){
+    m_write_type=WRITE_TYPE_RAW_WITH_PARTICLE;
+  }else if(savetype=="VTI"){
+    m_write_type=WRITE_TYPE_VTI;
+  }else if(savetype=="VTU"){
+    m_write_type=WRITE_TYPE_VTU;
+  /****fluid contents: end****/
   } else {
     cerr
       << "AFieldMaster: unknown output file format '"
@@ -99,9 +111,16 @@
   case WRITE_TYPE_RAW : suffix=".dat"; break;
   case WRITE_TYPE_RAW_WITH_ID : suffix=".dat"; break;
   case WRITE_TYPE_RAW_WITH_POS_ID : suffix=".dat"; break;
+  /****fluid contents: begin****/
+  case WRITE_TYPE_RAW_WITH_POS: suffix=".dat"; break; 
+  case WRITE_TYPE_RAW_WITH_ID_POS : suffix=".dat"; break;
+  case WRITE_TYPE_RAW_WITH_PARTICLE : suffix=".dat"; break; 
+  case WRITE_TYPE_VTI : suffix=".vti"; break;
+  case WRITE_TYPE_VTU : suffix=".vtu"; break;
+  /****fluid contents: end****/
   default : cerr << "AFieldMaster: wrong m_write_type in makeFilename" << endl; // can't happen
   }
-  fn << m_file_name << "." << m_save_count << suffix; 
+  fn << m_file_name << "." << m_save_count*m_dt << suffix;
   m_save_count++;
 
   return fn.str();
@@ -149,6 +168,13 @@
   case WRITE_TYPE_RAW2 : writeAsRAW2(); break;
   case WRITE_TYPE_RAW_WITH_ID : writeAsRawWithID(); break;
   case WRITE_TYPE_RAW_WITH_POS_ID : writeAsRawWithPosID(); break;
-  default : cerr << "AFieldMaster: wrong m_write_type in write" << endl; 
+  /***fluid contents:begin***/
+  case WRITE_TYPE_RAW_WITH_POS : writeAsRawWithPos(); break;
+  case WRITE_TYPE_RAW_WITH_ID_POS : writeAsRawWithIDPos(); break;
+  case WRITE_TYPE_RAW_WITH_PARTICLE : writeAsRawWithParticle(); break;
+  case WRITE_TYPE_VTI : writeAsVTI(); break;
+  case WRITE_TYPE_VTU : writeAsVTU(); break;
+  /***fluid contents:end***/
+  default : cerr << "AFieldMaster: wrong m_write_type in write" << endl;
   }
 }

=== modified file 'Fields/FieldMaster.h'
--- Fields/FieldMaster.h	2017-01-04 08:26:43 +0000
+++ Fields/FieldMaster.h	2019-07-28 14:55:57 +0000
@@ -40,7 +40,7 @@
 {
  private:
   static int s_field_count;
-  
+
  protected:
   TML_Comm *m_comm;
   string m_field_name;
@@ -64,6 +64,13 @@
   virtual void writeAsRAW(){console.Error()<<"writeAsRAW NOT IMPLEMENTED\n";};
   virtual void writeAsRawWithID(){console.Error()<<"writeAsRawWithID NOT IMPLEMENTED\n";};
   virtual void writeAsRawWithPosID(){console.Error()<<"writeAsRawWithPosID NOT IMPLEMENTED\n";};
+  /****fluid contents: begin***/
+  virtual void writeAsRawWithPos(){console.Error()<<"writeAsRawWithPosID NOT IMPLEMENTED\n";}; 
+  virtual void writeAsRawWithIDPos(){console.Error()<<"writeAsRawWithPosID NOT IMPLEMENTED\n";}; 
+  virtual void writeAsRawWithParticle(){console.Error()<<"writeAsRawWithPosID NOT IMPLEMENTED\n";}; 
+  virtual void writeAsVTI(){console.Error()<<"writeAsVTI NOT IMPLEMENTED\n";};
+  virtual void writeAsVTU(){console.Error()<<"writeAsVTU NOT IMPLEMENTED\n";};
+  /****fluid contents: end***/
 
  public:
   AFieldMaster(TML_Comm*,const string&,const string&,const string&,int,int,int);

=== added file 'Fields/FluidFieldMaster.cpp'
--- Fields/FluidFieldMaster.cpp	1970-01-01 00:00:00 +0000
+++ Fields/FluidFieldMaster.cpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,566 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+// --- TML includes ---
+#include "tml/comm/comm.h"
+
+//--- project includes ---
+#include "FieldMaster.h"
+#include "Foundation/vec3.h"
+#include "field_const.h"
+#include "FluidFieldMaster.h"
+
+// --- TML includes ---
+#include "tml/comm/comm.h"
+
+//--- IO includes ---
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+
+//--- STL inculdes ---
+#include <string>
+#include <map>
+
+
+using std::map;
+using std::multimap;
+using std::cout;
+using std::endl;
+using std::ofstream;
+using std::string;
+
+//==== SCALAR PFM ===
+
+/*!
+  Constructor. Set up the Master and broadcast parameters to the slaves.
+
+  \param comm the communicator
+  \param fieldname the name of the field to be saved
+  \param filename the name of the file to be saved into or the base for the generation of the filenames if the saving format requires multiple files
+  \param savetype the way to save data, currently supported are DX,SUM
+  \param t0 the first timestep to be saved
+  \param tend the last timestep to be saved
+  \param dt the interval between timesteps to be saving
+*/
+ScalarFluidFieldMaster::ScalarFluidFieldMaster(TML_Comm* comm,const string& fieldname,const string& filename,const string& savetype,int t0,int tend,int dt)
+  :AFieldMaster(comm,fieldname,filename,savetype,t0,tend,dt)
+{
+  m_comm->broadcast_cont(fieldname);
+  m_comm->broadcast(m_id);
+}
+
+
+void ScalarFluidFieldMaster::collect()
+{
+  // send field id to slave
+  m_comm->broadcast(m_id);
+
+  if((m_write_type==WRITE_TYPE_SUM)||(m_write_type==WRITE_TYPE_MAX)) {
+    collectSum();
+  } else {
+    collectFull();
+  }
+}
+
+/*!
+  collect the full set of data, i.e. value, and posn of fluid cells
+*/
+void ScalarFluidFieldMaster::collectFull()
+{
+  multimap<int,pair<Vec3,double> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=1;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<Vec3,double> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_save_vector.push_back(iter->second);
+  }
+}
+
+
+/*!
+  collect sum of values only
+*/
+void ScalarFluidFieldMaster::collectSum()
+{
+
+  multimap<int,double> temp_mm;
+
+  // send type of collect to slaves
+  m_comm->broadcast(m_write_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+
+  // collate receive data from multimap into single map
+  for(multimap<int,double>::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_sum_vector.push_back(iter->second);
+  }
+}
+
+
+/*!
+  sum data and write them out into a single continuous file
+
+  \warning n
+*/
+void ScalarFluidFieldMaster::writeAsSUM()
+{
+
+  // sum data
+  double sum_data=0.0;
+  for(vector<double>::iterator iter=m_sum_vector.begin();
+      iter!=m_sum_vector.end();
+      iter++){
+    sum_data+=*iter;
+  }
+  // open file for appending
+  ofstream out_file(m_file_name.c_str(),ios::app);
+  // write data
+  out_file << sum_data << endl;
+  // close file
+  out_file.close();
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+  m_sum_vector.erase(m_sum_vector.begin(),m_sum_vector.end());
+}
+
+/*!
+  find max datum and write it out into a single continuous file
+*/
+void ScalarFluidFieldMaster::writeAsMAX()
+{
+
+  // get max data
+  double max_data=*m_sum_vector.begin();
+  for(vector<double>::iterator iter=m_sum_vector.begin();
+      iter!=m_sum_vector.end();
+      iter++){
+    max_data=(*iter > max_data) ? *iter : max_data;
+  }
+  // open file for appending
+  ofstream out_file(m_file_name.c_str(),ios::app);
+  // write data
+  out_file << max_data << endl;
+  // close file
+  out_file.close();
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+  m_sum_vector.erase(m_sum_vector.begin(),m_sum_vector.end());
+}
+
+/*!
+  write data as a raw series of values, one row of values per time step,
+  all timesteps into the same file
+*/
+void ScalarFluidFieldMaster::writeAsRAW_SERIES()
+{
+
+  // open file
+  ofstream out_file(m_file_name.c_str(),ios::app);
+
+  // write data
+  for(vector<pair<Vec3,double> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    out_file << iter->first  << " " << iter->second << endl;
+  }
+
+  out_file << endl;
+  out_file.close();
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+  m_sum_vector.erase(m_sum_vector.begin(),m_sum_vector.end());
+}
+
+/*!
+  write data as raw position, value, one time step per file
+*/
+
+void ScalarFluidFieldMaster::writeAsRAW()
+{
+
+  //generate filename
+  string fn=makeFilename();
+
+  // open file
+  ofstream out_file(fn.c_str());
+
+  // write data
+
+  for(vector<pair<Vec3,double> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    out_file << iter->first  << " " << iter->second << endl;
+  }
+
+  // close file
+  out_file.close();
+
+  // clean up
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+  m_sum_vector.erase(m_sum_vector.begin(),m_sum_vector.end());
+}
+
+
+bool sortOnZ(const pair<Vec3,double> a, const pair<Vec3,double> b) {
+  Vec3 a_index=a.first;
+  Vec3 b_index=b.first;
+  return a_index.Z() < b_index.Z();
+}
+
+bool sortOnY(const pair<Vec3,double> a, const pair<Vec3,double> b) {
+  Vec3 a_index=a.first;
+  Vec3 b_index=b.first;
+  return a_index.Y() < b_index.Y();
+}
+
+bool sortOnX(const pair<Vec3,double> a, const pair<Vec3,double> b) {
+  Vec3 a_index=a.first;
+  Vec3 b_index=b.first;
+  return a_index.X() < b_index.X();
+}
+
+vector<pair<Vec3,double> > sortVector(vector<pair<Vec3,double> > unsorted) 
+{
+  vector<pair<Vec3,double> > sorted = unsorted;
+  std::stable_sort(sorted.begin(), sorted.end(), sortOnX);
+  std::stable_sort(sorted.begin(), sorted.end(), sortOnY);
+  std::stable_sort(sorted.begin(), sorted.end(), sortOnZ);
+  return sorted;
+}
+
+
+/*!
+  write data as VTI file for paraview, one time step per file
+*/
+
+void ScalarFluidFieldMaster::writeAsVTI()
+{
+
+  //generate filename
+  string fn=makeFilename();
+
+  // open file
+  ofstream out_file(fn.c_str());
+
+  vector<pair<Vec3,double> > sorted;
+  sorted=sortVector(m_save_vector);
+  double minX,minY,minZ,maxX,maxY,maxZ;
+  vector<pair<Vec3,double> >::iterator iter_begin=sorted.begin();
+  Vec3 pos_begin=iter_begin->first;
+  minX=pos_begin.X();minY=pos_begin.Y();minZ=pos_begin.Z();
+  vector<pair<Vec3,double> >::iterator iter_end=sorted.end()-1;
+  Vec3 pos_end=iter_end->first;
+  maxX=pos_end.X();maxY=pos_end.Y();maxZ=pos_end.Z();
+
+  int z_count=0;
+  for(vector<pair<Vec3,double> >::iterator iter0=sorted.begin();iter0!=sorted.end();iter0++){
+    Vec3 pos=iter0->first;
+    if(pos.Z()==minZ) {z_count++;} 
+    else {break;};
+  }
+  int nz=sorted.size()/z_count; 
+
+  int y_count=0;
+  for(vector<pair<Vec3,double> >::iterator iter0=sorted.begin();iter0!=sorted.end();iter0++){
+    Vec3 pos=iter0->first;
+    if(pos.Y()==minY) {y_count++;} 
+    else {break;};
+  }
+  int ny=z_count/y_count;
+
+  int x_count=0;
+  for(vector<pair<Vec3,double> >::iterator iter0=sorted.begin();iter0!=sorted.end();iter0++){
+    Vec3 pos=iter0->first;
+    if(pos.X()==minX) {x_count++;} 
+    else {break;};
+  }
+  int nx=y_count/x_count;
+
+  double xside,yside,zside;
+  if(nx==1){
+    xside=maxX-minX;
+  }else{
+    xside=(maxX-minX)/double(nx-1);
+  }
+
+  if(ny==1){
+    yside=maxY-minY;
+  }else{
+    yside=(maxY-minY)/double(ny-1);
+  }
+  if(nz==1){
+    zside=maxZ-minZ;
+  }else{
+    zside=(maxZ-minZ)/double(nz-1);
+  }
+
+  out_file<<"<VTKFile type=\"ImageData\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">"<< endl;
+  out_file<<"<ImageData WholeExtent=\"0 "<<nx<<" 0 "<<ny<<" 0 " <<nz<<"\" Origin=\""<<minX-0.5*xside<<" "<<minY-0.5*yside<<" "<<minZ-0.5*zside<<"\" Spacing=\""\
+          <<xside<<" "<<yside<<" "<<zside<<"\">" << endl;
+  out_file << "<Piece Extent=\"0 "<<nx<<" 0 "<<ny<<" 0 "<<nz<<"\">" << endl;
+  out_file << "<CellData>" << endl;
+  out_file << "<DataArray type=\"Float32\" Name=\"fluid scalar data\" format=\"ascii\">" << endl;
+
+  vector<pair<Vec3,double> >::iterator iter=sorted.begin();
+  for (int z=0;z<nz;z++){
+    for (int y=0;y<ny;y++){
+      for (int x=0;x<nx;x++){
+        out_file << iter->second << " ";
+        iter++;
+      };
+      out_file<<endl;
+    };
+    out_file<<endl;  
+  };  
+
+  out_file << "</DataArray>" << endl <<endl;
+  out_file << "</CellData>" << endl;
+  out_file << "</Piece>" << endl;
+  out_file << "</ImageData>" << endl;
+  out_file << "</VTKFile>" << endl;
+
+  // close file
+  out_file.close();
+
+  // clean up
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+  m_sum_vector.erase(m_sum_vector.begin(),m_sum_vector.end());
+}
+
+
+// === VECTOR PFM ===
+/*!
+  Constructor. Set up the Master and broadcast parameters to the slaves.
+
+  \param comm the communicator
+  \param fieldname the name of the field to be saved
+  \param filename the name of the file to be saved into or the base for the generation of the filenames if the saving format requires multiple files
+  \param savetype the way to save data, currently supported are DX,SUM
+  \param t0 the first timestep to be saved
+  \param tend the last timestep to be saved
+  \param dt the interval between timesteps to be saving
+*/
+VectorFluidFieldMaster::VectorFluidFieldMaster(TML_Comm* comm,const string& fieldname,const string& filename,const string& savetype,int t0,int tend,int dt)
+  :AFieldMaster(comm,fieldname,filename,savetype,t0,tend,dt)
+{
+  m_comm->broadcast_cont(fieldname);
+  m_comm->broadcast(m_id);
+}
+
+void VectorFluidFieldMaster::collect()
+{
+  multimap<int,pair<Vec3,Vec3> > temp_mm;
+
+  // send field id to slaves
+  m_comm->broadcast(m_id);
+
+  //receive data from fields
+  m_comm->gather(temp_mm);
+
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<Vec3,Vec3> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_save_vector.push_back(iter->second);
+  }
+}
+
+
+/*!
+  sum data and write them out into a single continuous file
+
+  \warning untested
+*/
+void VectorFluidFieldMaster::writeAsSUM()
+{
+  Vec3 sum_data=Vec3(0.0,0.0,0.0);
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    sum_data+=iter->second;
+  }
+  // open file for appending
+  ofstream out_file(m_file_name.c_str(),ios::app);
+  // write data
+  out_file << sum_data << endl;
+  // close file
+  out_file.close();
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+}
+
+/*!
+  get max datum and write it out into a single continuous file
+
+  \warning untested
+*/
+void VectorFluidFieldMaster::writeAsMAX()
+{
+  Vec3 max_data=(m_save_vector.begin())->second;
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    max_data=(iter->second.norm2() > max_data.norm2()) ? iter->second : max_data;
+  }
+  // open file for appending
+  ofstream out_file(m_file_name.c_str(),ios::app);
+  // write data
+  out_file << max_data << endl;
+  // close file
+  out_file.close();
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+}
+
+
+/*!
+  write data as a raw series of values, one row of values per time step,
+  all timesteps into the same file
+*/
+void VectorFluidFieldMaster::writeAsRAW_SERIES()
+{
+  // open file
+  ofstream out_file(m_file_name.c_str(),ios::app);
+
+  // write data
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    out_file << iter->first  << " " << iter->second << endl;
+  }
+  out_file << endl;
+  out_file.close();
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+}
+
+/*!
+  write data as raw position, value pairs, one time step per file
+*/
+void VectorFluidFieldMaster::writeAsRAW()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // open file
+  ofstream out_file(fn.c_str());
+
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    out_file << iter->first  << " " << iter->second << endl;
+  }
+
+  // close file
+  out_file.close();
+
+  // clean up
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+}
+
+/*!
+  write data as VTU file for paraview, position, value, one time step per file
+*/
+void VectorFluidFieldMaster::writeAsVTU()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // open file
+  ofstream out_file(fn.c_str());
+
+  int n=0;
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    n++;
+  }
+  double xside,yside,zside;
+  Vec3 number0;
+  vector<pair<Vec3,Vec3> >::iterator iter0=m_save_vector.begin();
+  Vec3 pos0=iter0->first;
+
+  for(iter0=m_save_vector.begin();iter0!=m_save_vector.end();iter0++){
+    number0=iter0->first;
+    if(number0.X()!=pos0.X()){xside=fabs(number0.X()-pos0.X());break;};
+  }
+  for(iter0=m_save_vector.begin();iter0!=m_save_vector.end();iter0++){
+    number0=iter0->first;
+    if(number0.Y()!=pos0.Y()){yside=fabs(number0.Y()-pos0.Y());break;};
+  }
+  for(iter0=m_save_vector.begin();iter0!=m_save_vector.end();iter0++){
+    number0=iter0->first;
+    if(number0.Z()!=pos0.Z()){zside=fabs(number0.Z()-pos0.Z());break;};
+  }
+
+  double cellsize=(xside+yside+zside)/3.0;
+  out_file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">" << endl;
+  out_file << "<UnstructuredGrid>" << endl;
+  out_file << "<Piece NumberOfPoints=\"" << n << "\" NumberOfCells=\"" << 0 << "\">" << endl;
+
+  //write positions
+  out_file << "<Points>" << endl;
+  out_file << "<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">" << endl;
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    out_file << iter->first<< endl;
+  }
+  out_file << "</DataArray>" << endl;
+  out_file << "</Points>" << endl;
+
+  //write cell size
+  out_file << "<PointData Scalars=\"cellsize\">" << endl;
+  out_file << "<DataArray type=\"Float64\" Name=\"cellsize\" NumberOfComponents=\"1\" format=\"ascii\">" << endl;
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    out_file << cellsize << endl;
+  }
+  out_file << "</DataArray>" << endl;
+
+  //write vector value
+  out_file << "<DataArray type=\"Float64\" Name=\"vector value\" NumberOfComponents=\"3\" format=\"ascii\">" << endl;
+  for(vector<pair<Vec3,Vec3> >::iterator iter=m_save_vector.begin();
+      iter!=m_save_vector.end();
+      iter++){
+    out_file << iter->second << endl;
+  }
+
+  out_file << "</DataArray>" << endl;
+  out_file << "</PointData>" << endl;
+  out_file << "<Cells>" << endl;
+  out_file << "<DataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"connectivity\" format=\"ascii\">" << endl;
+  out_file << "</DataArray>" << endl;
+  out_file << "<DataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"offsets\" format=\"ascii\">" << endl;
+  out_file << "</DataArray>" << endl;
+  out_file << "<DataArray type=\"UInt8\" NumberOfComponents=\"1\" Name=\"types\" format=\"ascii\">" << endl;
+  out_file << "</DataArray>" << endl;
+  out_file << "</Cells>" << endl;
+  out_file << "</Piece>" << endl;
+  out_file << "</UnstructuredGrid>" << endl;
+  out_file << "</VTKFile>" << endl;
+
+  // close file
+  out_file.close();
+
+  // clean up
+  m_save_vector.erase(m_save_vector.begin(),m_save_vector.end());
+}
+
+

=== added file 'Fields/FluidFieldMaster.h'
--- Fields/FluidFieldMaster.h	1970-01-01 00:00:00 +0000
+++ Fields/FluidFieldMaster.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,106 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+#ifndef __FLUIDFIELDMASTER_H
+#define __FLUIDFIELDMASTER_H
+
+//--- project includes ---
+#include "FieldMaster.h"
+#include "Foundation/vec3.h"
+#include "field_const.h"
+
+// --- TML includes ---
+#include "tml/comm/comm.h"
+
+//--- IO includes ---
+#include <iostream>
+#include <fstream>
+
+//--- STL inculdes ---
+#include <string>
+#include <map>
+
+
+using std::map;
+using std::multimap;
+using std::cout;
+using std::endl;
+using std::ofstream;
+using std::string;
+
+
+class TML_Comm;
+
+bool sortOnZ(const pair<Vec3,double>, const pair<Vec3,double>); 
+bool sortOnY(const pair<Vec3,double>, const pair<Vec3,double>); 
+bool sortOnX(const pair<Vec3,double>, const pair<Vec3,double>);
+vector<pair<Vec3,double> > sortVector(vector<pair<Vec3,double> >);
+
+/*!
+  \class ScalarFluidFieldMaster
+  \brief Class for master part of a scalar field which is defined on all fluid cells
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+class ScalarFluidFieldMaster : public AFieldMaster
+{
+ protected:
+  vector<pair<Vec3,double> > m_save_vector;
+  vector<double> m_sum_vector;
+
+  void collectFull();
+  void collectSum();
+  virtual void writeAsSUM();
+  virtual void writeAsMAX();
+  virtual void writeAsRAW_SERIES();
+  virtual void writeAsRAW(); 
+  virtual void writeAsVTI();
+
+ public:
+  ScalarFluidFieldMaster(TML_Comm*,const string&,const string&,const string&,int,int,int);
+  virtual ~ScalarFluidFieldMaster(){};
+
+  virtual void collect();
+
+};
+
+/*!
+  \class VectorFluidFieldMaster
+  \brief Class for master part of a vector field which is defined on all fluid cells
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+class VectorFluidFieldMaster : public AFieldMaster
+{
+ protected:
+  vector<pair<Vec3,Vec3> > m_save_vector;
+
+  virtual void writeAsSUM();
+  virtual void writeAsMAX();
+  virtual void writeAsRAW_SERIES();
+  virtual void writeAsRAW();
+  virtual void writeAsVTU();
+
+ public:
+  VectorFluidFieldMaster(TML_Comm*,const string&,const string&,const string&,int,int,int);
+
+  virtual ~VectorFluidFieldMaster(){};
+
+  void collect();
+};
+
+
+#endif //__FLUIDFIELDMASTER_H

=== added file 'Fields/FluidInteractionFieldMaster.cpp'
--- Fields/FluidInteractionFieldMaster.cpp	1970-01-01 00:00:00 +0000
+++ Fields/FluidInteractionFieldMaster.cpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,695 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+// --- TML includes ---
+#include "tml/comm/comm.h"
+
+// --- Project includes ---
+#include "FluidInteractionFieldMaster.h"
+#include "field_const.h"
+#include "console.h"
+#include "FieldMaster.h"
+#include "Foundation/vec3.h"
+#include "field_const.h"
+
+
+//--- IO includes ---
+#include <iostream>
+#include <fstream>
+
+//--- STL inculdes ---
+#include <string>
+#include <map>
+
+using std::cout;
+using std::endl;
+using std::ofstream;
+using std::string;
+using std::multimap;
+
+
+//==== SCALAR PFM ===
+/*!
+  Constructor. Setup master and broadcast parameters to slaves
+
+  \param comm the communicator
+  \param fieldname the name of the field to be saved
+  \param filename the name of the file to be saved into or the base for the generation of the filenames if the saving format requires multiple files
+  \param savetype the way to save data, currently supported are DX,SUM
+  \param t0 the first timestep to be saved
+  \param tend the last timestep to be saved
+  \param dt the interval between timesteps to be saving
+  \param checked choose between "full" and "checked" fields
+*/
+ScalarFluidInteractionFieldMaster::ScalarFluidInteractionFieldMaster(TML_Comm* comm,const string& fieldname,const string& filename,const string& savetype,int t0,int tend,int dt)
+  :AFieldMaster(comm,fieldname,filename,savetype,t0,tend,dt)
+{
+  m_comm->broadcast_cont(fieldname);
+  m_comm->broadcast(m_id);
+}
+
+
+void ScalarFluidInteractionFieldMaster::collect()
+{
+  // send field id to slave
+  m_comm->broadcast(m_id);
+
+  switch(m_write_type){
+  case WRITE_TYPE_SUM: collectSum(); break;
+  case WRITE_TYPE_MAX: collectMax(); break;
+  case WRITE_TYPE_RAW_WITH_POS: collectFullwithPos(); break;
+  case WRITE_TYPE_RAW_WITH_ID: collectFullwithID(); break;
+  case WRITE_TYPE_RAW_WITH_ID_POS: collectFullwithIDPos(); break;
+  case WRITE_TYPE_RAW_WITH_PARTICLE: collectFullwithParticle(); break;
+
+  default: collectFullwithParticle();
+  }
+}
+
+
+/*!
+  collect full data set, particle position
+*/
+void ScalarFluidInteractionFieldMaster::collectFullwithPos()
+{
+  multimap<int,pair<Vec3,double> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_POS;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<Vec3,double> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_pos.push_back(iter->second);
+  }
+}
+
+/*!
+  collect full data set, particle id
+*/
+void ScalarFluidInteractionFieldMaster::collectFullwithID()
+{
+  multimap<int,pair<int,double> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_ID;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<int,double> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_id.push_back(iter->second);
+  }
+}
+
+/*!
+  collect full data set, particle id and position
+*/
+void ScalarFluidInteractionFieldMaster::collectFullwithIDPos()
+{
+  multimap<int,pair<pair<int,Vec3>,double> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_ID_POS;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<pair<int,Vec3>,double> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_id_pos.push_back(iter->second);
+  }
+}
+
+/*!
+  collect full data set, particle id, position and radius
+*/
+void ScalarFluidInteractionFieldMaster::collectFullwithParticle()
+{
+  multimap<int,pair<esys::lsm::triplet<int,Vec3,double>,double> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_PARTICLE;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<esys::lsm::triplet<int,Vec3,double>,double> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_particle.push_back(iter->second);
+  }
+}
+
+/*!
+  collect sum of values only
+*/
+void ScalarFluidInteractionFieldMaster::collectSum()
+{
+  multimap<int,double> temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_SUM;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+
+  // collate receive data from multimap into single map
+  for(multimap<int,double>::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_sum_vec.push_back(iter->second);
+  }
+}
+
+/*!
+  collect max of values only
+*/
+void ScalarFluidInteractionFieldMaster::collectMax()
+{
+  multimap<int,double> temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_MAX;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+
+  // collate receive data from multimap into single map
+  for(multimap<int,double>::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_max_vec.push_back(iter->second);
+  }
+}
+
+
+/*!
+  write data as pos,value groups
+*/
+void ScalarFluidInteractionFieldMaster::writeAsRawWithPos()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // write data
+  ofstream out_file(fn.c_str());
+  // check if file is sucessfully opened
+  if(!out_file){
+    console.Error() << "can not open output file " << fn << "\n";
+  } else {
+    int count=0;
+    console.XDebug() << m_data_with_pos.size() << " vectors to be written\n";
+    for(vector<pair<Vec3,double> >::iterator iter=m_data_with_pos.begin();
+        iter!=m_data_with_pos.end();
+        iter++) {
+      out_file << iter->first << " " << iter->second << endl;
+      count++;
+      // write debug message every 10k writes
+      if((count % 10000)==0){
+        console.XDebug() << count << " vectors written\n";
+      }
+    }
+    console.XDebug() << "finished writing " << count << " vectors \n";
+    // close file
+    out_file.close();
+  }
+  //clean up
+  m_data_with_pos.clear();
+}
+
+/*!
+  write data as id,value groups
+*/
+void ScalarFluidInteractionFieldMaster::writeAsRawWithID()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // write data
+  ofstream out_file(fn.c_str());
+  // check if file is sucessfully opened
+  if(!out_file){
+    console.Error() << "can not open output file " << fn << "\n";
+  } else {
+    int count=0;
+    console.XDebug() << m_data_with_id.size() << " vectors to be written\n";
+    for(vector<pair<int,double> >::iterator iter=m_data_with_id.begin();
+        iter!=m_data_with_id.end();
+        iter++) {
+      out_file << iter->first << " " << iter->second << endl;
+      count++;
+      // write debug message every 10k writes
+      if((count % 10000)==0){
+        console.XDebug() << count << " vectors written\n";
+      }
+    }
+    console.XDebug() << "finished writing " << count << " vectors \n";
+    // close file
+    out_file.close();
+  }
+  //clean up
+  m_data_with_id.clear();
+}
+
+/*!
+  write data as id,pos,value groups
+*/
+void ScalarFluidInteractionFieldMaster::writeAsRawWithIDPos()
+{
+  cerr << "ScalarFluidInteractionFieldMaster::writeAsRawWithIDPos not implemented" << endl;
+}
+
+/*!
+  write data as id,pos,radius,value groups
+*/
+void ScalarFluidInteractionFieldMaster::writeAsRawWithParticle()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // write data
+  ofstream out_file(fn.c_str());
+  // check if file is sucessfully opened
+  if(!out_file){
+    console.Error() << "can not open output file " << fn << "\n";
+  } else {
+    int count=0;
+    console.XDebug() << m_data_with_particle.size() << " vectors to be written\n";
+    for(vector<pair<esys::lsm::triplet<int,Vec3,double>,double> >::iterator iter=m_data_with_particle.begin();
+        iter!=m_data_with_particle.end();
+        iter++) {
+      out_file
+        << (iter->first).get<0>()
+        << " "
+        << (iter->first).get<1>()
+        << " "
+        << (iter->first).get<2>()
+	<< " "
+        << iter->second
+        << endl;
+      count++;
+      // write debug message every 10k writes
+      if((count % 10000)==0){
+        console.XDebug() << count << " vectors written\n";
+      }
+    }
+    console.XDebug() << "finished writing " << count << " vectors \n";
+    // close file
+    out_file.close();
+  }
+  //clean up
+  m_data_with_particle.clear();
+}
+
+
+/*!
+  sum data and write them out into a single continuous file
+*/
+void ScalarFluidInteractionFieldMaster::writeAsSUM()
+{
+  // sum data
+  double sum_data=0.0;
+  for(vector<double>::iterator iter=m_sum_vec.begin();
+      iter!=m_sum_vec.end();
+      iter++){
+    sum_data+=*iter;
+  }
+  // open file for appending
+  ofstream out_file(m_file_name.c_str(),ios::app);
+
+  // set output precision to 10 significant digits, scientific format
+  out_file.setf(ios_base::scientific, ios_base::floatfield);
+  out_file.precision(10);
+
+  // write data
+  out_file << sum_data << endl;
+
+  // close file
+  out_file.close();
+  //clean up
+  m_sum_vec.erase(m_sum_vec.begin(),m_sum_vec.end());
+
+}
+
+/*!
+  get the maximum of the data and write it out into a single continuous file
+*/
+void ScalarFluidInteractionFieldMaster::writeAsMAX()
+{
+  // sum data
+  double max_data=*(m_max_vec.begin());
+  for(vector<double>::iterator iter=m_max_vec.begin();
+      iter!=m_max_vec.end();
+      iter++){
+    max_data=(*iter > max_data) ? *iter : max_data;
+  }
+  // open file for appending
+  ofstream out_file(m_file_name.c_str(),ios::app);
+  // write data
+  out_file << max_data << endl;
+
+  // close file
+  out_file.close();
+  //clean up
+  m_max_vec.erase(m_max_vec.begin(),m_max_vec.end());
+
+}
+
+/*!
+  write data as a raw series of values, one row of values per time step,
+  all timesteps into the same file
+*/
+void ScalarFluidInteractionFieldMaster::writeAsRAW_SERIES()
+{
+  cerr << "ScalarFluidInteractionFieldMaster::writeAsRAW_SERIES not implemented" << endl;
+}
+
+
+
+
+
+// === VECTOR PFM ===
+/*!
+  Constructor. Setup master and broadcast parameters to slaves
+
+  \param comm the communicator
+  \param fieldname the name of the field to be saved
+  \param filename the name of the file to be saved into or the base for the generation of the filenames if the saving format requires multiple files
+  \param savetype the way to save data, currently supported are DX,SUM
+  \param t0 the first timestep to be saved
+  \param tend the last timestep to be saved
+  \param dt the interval between timesteps to be saving
+  \param checked choose between "full" and "checked" fields
+*/
+VectorFluidInteractionFieldMaster::VectorFluidInteractionFieldMaster(TML_Comm* comm,const string& fieldname,const string& filename,const string& savetype,int t0,int tend,int dt)
+  :AFieldMaster(comm,fieldname,filename,savetype,t0,tend,dt)
+{
+  m_comm->broadcast_cont(fieldname);
+  m_comm->broadcast(m_id);
+}
+
+
+void VectorFluidInteractionFieldMaster::collect()
+{
+  // send field id to slave
+  m_comm->broadcast(m_id);
+
+  switch(m_write_type){
+  case WRITE_TYPE_SUM: collectSum(); break;
+  //case WRITE_TYPE_MAX: collectMax(); break;
+  case WRITE_TYPE_RAW_WITH_POS: collectFullwithPos(); break;
+  case WRITE_TYPE_RAW_WITH_ID: collectFullwithID(); break;
+  case WRITE_TYPE_RAW_WITH_ID_POS: collectFullwithIDPos(); break;
+  case WRITE_TYPE_RAW_WITH_PARTICLE: collectFullwithParticle(); break;
+
+  default: collectFullwithParticle();
+  }
+}
+
+
+/*!
+  collect full data set, particle position
+*/
+void VectorFluidInteractionFieldMaster::collectFullwithPos()
+{
+  multimap<int,pair<Vec3,Vec3 > > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_POS;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<Vec3,Vec3> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_pos.push_back(iter->second);
+  }
+}
+
+/*!
+  collect full data set, particle id
+*/
+void VectorFluidInteractionFieldMaster::collectFullwithID()
+{
+  multimap<int,pair<int,Vec3> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_ID;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<int,Vec3> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_id.push_back(iter->second);
+  }
+}
+
+/*!
+  collect full data set, particle id and position
+*/
+void VectorFluidInteractionFieldMaster::collectFullwithIDPos()
+{
+  multimap<int,pair<pair<int,Vec3>,Vec3> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_ID_POS;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<pair<int,Vec3>,Vec3> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_id_pos.push_back(iter->second);
+  }
+}
+
+/*!
+  collect full data set, particle id, position and radius
+*/
+void VectorFluidInteractionFieldMaster::collectFullwithParticle()
+{
+  multimap<int,pair<esys::lsm::triplet<int,Vec3,double>,Vec3> > temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_FULL_WITH_PARTICLE;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+  // collate receive data from multimap into single map
+  for(multimap<int,pair<esys::lsm::triplet<int,Vec3,double>,Vec3> >::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_data_with_particle.push_back(iter->second);
+  }
+}
+
+/*!
+  collect sum of values only
+*/
+void VectorFluidInteractionFieldMaster::collectSum()
+{
+  multimap<int,Vec3> temp_mm;
+
+  // send type of collect to slaves
+  int coll_type=COLL_TYPE_SUM;
+  m_comm->broadcast(coll_type);
+
+  // get data from slaves
+  m_comm->gather(temp_mm);
+
+  // collate receive data from multimap into single map
+  for(multimap<int,Vec3>::iterator iter=temp_mm.begin();
+      iter!=temp_mm.end();
+      iter++){
+    m_sum_vec.push_back(iter->second);
+  }
+}
+
+
+/*!
+  write data as pos,value groups
+*/
+void VectorFluidInteractionFieldMaster::writeAsRawWithPos()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // write data
+  ofstream out_file(fn.c_str());
+  // check if file is sucessfully opened
+  if(!out_file){
+    console.Error() << "can not open output file " << fn << "\n";
+  } else {
+    int count=0;
+    console.XDebug() << m_data_with_pos.size() << " vectors to be written\n";
+    for(vector<pair<Vec3,Vec3> >::iterator iter=m_data_with_pos.begin();
+        iter!=m_data_with_pos.end();
+        iter++) {
+      out_file << iter->first << " " << iter->second << endl;
+      count++;
+      // write debug message every 10k writes
+      if((count % 10000)==0){
+        console.XDebug() << count << " vectors written\n";
+      }
+    }
+    console.XDebug() << "finished writing " << count << " vectors \n";
+    // close file
+    out_file.close();
+  }
+  //clean up
+  m_data_with_pos.clear();
+}
+
+/*!
+  write data as id,value groups
+*/
+void VectorFluidInteractionFieldMaster::writeAsRawWithID()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // write data
+  ofstream out_file(fn.c_str());
+  // check if file is sucessfully opened
+  if(!out_file){
+    console.Error() << "can not open output file " << fn << "\n";
+  } else {
+    int count=0;
+    console.XDebug() << m_data_with_id.size() << " vectors to be written\n";
+    for(vector<pair<int,Vec3> >::iterator iter=m_data_with_id.begin();
+        iter!=m_data_with_id.end();
+        iter++) {
+      out_file << iter->first << " " << iter->second << endl;
+      count++;
+      // write debug message every 10k writes
+      if((count % 10000)==0){
+        console.XDebug() << count << " vectors written\n";
+      }
+    }
+    console.XDebug() << "finished writing " << count << " vectors \n";
+    // close file
+    out_file.close();
+  }
+  //clean up
+  m_data_with_id.clear();
+}
+
+/*!
+  write data as id,pos,value groups
+*/
+void VectorFluidInteractionFieldMaster::writeAsRawWithIDPos()
+{
+  cerr << "VectorFluidInteractionFieldMaster::writeAsRawWithIDPos not implemented" << endl;
+}
+
+/*!
+  write data as id,pos,radius,value groups
+*/
+void VectorFluidInteractionFieldMaster::writeAsRawWithParticle()
+{
+  //generate filename
+  string fn=makeFilename();
+
+  // write data
+  ofstream out_file(fn.c_str());
+  // check if file is sucessfully opened
+  if(!out_file){
+    console.Error() << "can not open output file " << fn << "\n";
+  } else {
+    int count=0;
+    console.XDebug() << m_data_with_particle.size() << " vectors to be written\n";
+    for(vector<pair<esys::lsm::triplet<int,Vec3,double>,Vec3> >::iterator iter=m_data_with_particle.begin();
+        iter!=m_data_with_particle.end();
+        iter++) {
+      out_file
+        << (iter->first).get<0>()
+        << " "
+        << (iter->first).get<1>()
+        << " "
+        << (iter->first).get<2>()
+	<< " "
+        << iter->second
+        << endl;
+      count++;
+      // write debug message every 10k writes
+      if((count % 10000)==0){
+        console.XDebug() << count << " vectors written\n";
+      }
+    }
+    console.XDebug() << "finished writing " << count << " vectors \n";
+    // close file
+    out_file.close();
+  }
+  //clean up
+  m_data_with_particle.clear();
+}
+
+
+/*!
+  sum data and write them out into a single continuous file
+*/
+void VectorFluidInteractionFieldMaster::writeAsSUM()
+{
+  // sum data
+  Vec3 sum_data=Vec3(0.0,0.0,0.0);
+  for(vector<Vec3>::iterator iter=m_sum_vec.begin();
+      iter!=m_sum_vec.end();
+      iter++){
+    sum_data+=*iter;
+  }
+  // open file for appending
+  ofstream out_file(m_file_name.c_str(),ios::app);
+
+  // set output precision to 10 significant digits, scientific format
+  out_file.setf(ios_base::scientific, ios_base::floatfield);
+  out_file.precision(10);
+
+  // write data
+  out_file << sum_data << endl;
+
+  // close file
+  out_file.close();
+  //clean up
+  m_sum_vec.erase(m_sum_vec.begin(),m_sum_vec.end());
+
+}
+
+/*!
+  write data as a raw series of values, one row of values per time step,
+  all timesteps into the same file
+*/
+void VectorFluidInteractionFieldMaster::writeAsRAW_SERIES()
+{
+  cerr << "VectorFluidInteractionFieldMaster::writeAsRAW_SERIES not implemented" << endl;
+}
+

=== added file 'Fields/FluidInteractionFieldMaster.h'
--- Fields/FluidInteractionFieldMaster.h	1970-01-01 00:00:00 +0000
+++ Fields/FluidInteractionFieldMaster.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,124 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+#ifndef __FLUIDINTERACTIONFIELDMASTER_H
+#define __FLUIDINTERACTIONFIELDMASTER_H
+
+//--- project includes ---
+#include "FieldMaster.h"
+#include "Foundation/vec3.h"
+#include "Foundation/quintuple.h"
+#include "Foundation/triplet.h"
+#include "field_const.h"
+
+// --- TML includes ---
+#include "tml/comm/comm.h"
+
+//--- IO includes ---
+#include <iostream>
+#include <fstream>
+
+//--- STL includes ---
+#include <vector>
+#include <string>
+#include <map>
+
+using std::vector;
+using std::map;
+using std::multimap;
+using std::cout;
+using std::endl;
+using std::ofstream;
+using std::string;
+
+class TML_Comm;
+
+/*!
+  \class ScalarFluidInteractionFieldMaster
+  \brief Class for master part of a scalar field which is defined on all fluid interactions
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+class ScalarFluidInteractionFieldMaster : public AFieldMaster
+{
+ protected:
+  vector<pair<Vec3,double> >  m_data_with_pos; // vector of <position,value> pairs
+  vector<pair<int,double> >  m_data_with_id;
+  vector<pair<pair<int,Vec3>,double> >  m_data_with_id_pos;
+  vector<pair<esys::lsm::triplet<int,Vec3,double>,double> > m_data_with_particle;
+
+  vector<double> m_sum_vec;
+  vector<double> m_max_vec;
+
+  virtual void writeAsRawWithPos();
+  virtual void writeAsRawWithID();
+  virtual void writeAsRawWithIDPos();
+  virtual void writeAsRawWithParticle();
+  virtual void writeAsSUM();
+  virtual void writeAsMAX();
+  virtual void writeAsRAW_SERIES();
+
+  void collectFullwithPos();
+  void collectFullwithID();
+  void collectFullwithIDPos();
+  void collectFullwithParticle();
+  void collectSum();
+  void collectMax();
+
+ public:
+  ScalarFluidInteractionFieldMaster(TML_Comm*,const string&,const string&,const string&,int,int,int);
+  virtual ~ScalarFluidInteractionFieldMaster(){};
+
+  void collect();
+};
+
+
+/*!
+  \class VectorFluidInteractionFieldMaster
+  \brief Class for master part of a vector field which is defined on all fluid interactions
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+class VectorFluidInteractionFieldMaster : public AFieldMaster
+{
+ protected:
+  vector<pair<Vec3,Vec3> > m_data_with_pos; // vector of <position,value> pairs
+  vector<pair<int,Vec3> >  m_data_with_id;
+  vector<pair<pair<int,Vec3>,Vec3> > m_data_with_id_pos;
+  vector<pair<esys::lsm::triplet<int,Vec3,double>,Vec3> > m_data_with_particle;
+
+  vector<Vec3> m_sum_vec;
+
+  virtual void writeAsRawWithPos();
+  virtual void writeAsRawWithID();
+  virtual void writeAsRawWithIDPos();
+  virtual void writeAsRawWithParticle();
+  virtual void writeAsSUM();
+  virtual void writeAsRAW_SERIES();
+
+  void collectFullwithPos();
+  void collectFullwithID();
+  void collectFullwithIDPos();
+  void collectFullwithParticle();
+  void collectSum();
+
+ public:
+  VectorFluidInteractionFieldMaster(TML_Comm*,const string&,const string&,const string&,int,int,int);
+  virtual ~VectorFluidInteractionFieldMaster(){};
+
+  void collect();
+};
+#endif //__FlUIDINTERACTIONFIELDMASTER_H

=== modified file 'Fields/Makefile.am'
--- Fields/Makefile.am	2017-01-04 08:26:43 +0000
+++ Fields/Makefile.am	2019-07-28 14:55:57 +0000
@@ -61,7 +61,19 @@
 	VectorEdge2DFieldSlave.h \
 	VectorEdge2DFieldSlave.cpp \
 	MaxTrigger.cpp \
-	MaxTrigger.h
+	MaxTrigger.h \
+	FluidFieldMaster.h \
+	FluidFieldMaster.cpp \
+	ScalarFluidFieldSlave.hpp \
+	ScalarFluidFieldSlave.h \
+	VectorFluidFieldSlave.hpp \
+	VectorFluidFieldSlave.h \
+	FluidInteractionFieldMaster.h \
+	FluidInteractionFieldMaster.cpp \
+	ScalarFluidInteractionFieldSlave.hpp \
+	ScalarFluidInteractionFieldSlave.h \
+	VectorFluidInteractionFieldSlave.hpp \
+	VectorFluidInteractionFieldSlave.h
 
 libFields_la_LIBADD =                           \
   $(top_builddir)/Foundation/libFoundation.la   \

=== added file 'Fields/ScalarFluidFieldSlave.h'
--- Fields/ScalarFluidFieldSlave.h	1970-01-01 00:00:00 +0000
+++ Fields/ScalarFluidFieldSlave.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,51 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+#ifndef __SCALAR_FLUID_FIELD_SLAVE_H
+#define __SCALAR_FLUID_FIELD_SLAVE_H
+
+// -- project includes --
+#include "FieldSlave.h"
+#include "Model/FluidCell.h"
+
+template <class T> class ParallelParticleArray;
+class TML_Comm;
+
+/*!
+  \class ScalarFluidFieldSlave
+  \brief class for slave part of scalar field defined on the fluid cells
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+template <typename T>
+class ScalarFluidFieldSlave : public AFieldSlave
+{
+ private:
+  virtual void SendDataFull();
+  virtual void SendDataSum();
+  virtual void SendDataMax();
+
+ protected:
+  CFluidCell::ScalarFieldFunction m_rdf;
+  ParallelParticleArray<T>* m_ppa;
+
+ public:
+  ScalarFluidFieldSlave(TML_Comm*,ParallelParticleArray<T>*,CFluidCell::ScalarFieldFunction);
+
+  virtual void sendData();
+};
+
+#include "ScalarFluidFieldSlave.hpp"
+
+#endif //__SCALAR_FLUID_FIELD_SLAVE_H

=== added file 'Fields/ScalarFluidFieldSlave.hpp'
--- Fields/ScalarFluidFieldSlave.hpp	1970-01-01 00:00:00 +0000
+++ Fields/ScalarFluidFieldSlave.hpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,121 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+//-- STL includes --
+#include <vector>
+#include <utility>
+
+using std::vector;
+using std::pair;
+
+// -- IO includes --
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+/*!
+  constructor
+
+  \param comm the TML communicator used for sending the data back to the master
+  \param ppa a pointer to the particle array
+  \param rdf the fluid cell member function to access the data
+*/
+template <typename T>
+ScalarFluidFieldSlave<T>::ScalarFluidFieldSlave(TML_Comm* comm,ParallelParticleArray<T>* ppa,CFluidCell::ScalarFieldFunction rdf):AFieldSlave(comm)
+{
+  m_ppa=ppa;
+  m_rdf=rdf;
+}
+
+
+
+/*!
+  send full field data and position of the fluid cells
+*/
+template <typename T>
+void ScalarFluidFieldSlave<T>::SendDataFull()
+{
+  vector<pair<Vec3,double> > data_vec;
+
+  data_vec=m_ppa->forAllInnerCellsGet(m_rdf);
+
+  // send data to master
+  m_comm->send_gather(data_vec,0);
+}
+
+/*!
+  send sum only
+*/
+template <typename T>
+void ScalarFluidFieldSlave<T>::SendDataSum()
+{
+  vector<double> data_vec;
+
+  // get data from particles
+  data_vec=m_ppa->forAllInnerCellsGetSum(m_rdf);
+
+  // sum data
+  double sum=0.0;
+  for(vector<double>::iterator iter=data_vec.begin();
+      iter!=data_vec.end();
+      iter++){
+    sum+=*iter;
+  }
+
+  vector<double> sum_vec;
+  sum_vec.push_back(sum);
+  m_comm->send_gather(sum_vec,0);
+}
+
+
+/*!
+  send maximum only
+*/
+template <typename T>
+void ScalarFluidFieldSlave<T>::SendDataMax()
+{
+  vector<double> data_vec;
+
+  // get data from particles
+  data_vec=m_ppa->forAllInnerCellsGetSum(m_rdf);
+
+  // sum data
+  double max=*(data_vec.begin());
+  for(vector<double>::iterator iter=data_vec.begin();
+      iter!=data_vec.end();
+      iter++){
+    max=(*iter > max) ? *iter : max;
+  }
+
+  vector<double> max_vec;
+  max_vec.push_back(max);
+  m_comm->send_gather(max_vec,0);
+}
+
+
+/*!
+  send data back to master
+*/
+template <typename T>
+void ScalarFluidFieldSlave<T>::sendData()
+{
+  int coll_type;
+  m_comm->recv_broadcast(coll_type,0);
+
+  switch(coll_type){
+  case 1: SendDataFull();break;
+  case 2: SendDataSum();break;
+  case 3: SendDataMax();break;
+  default: std::cerr << "unknown collection type" << std::endl;
+  }
+}

=== added file 'Fields/ScalarFluidInteractionFieldSlave.h'
--- Fields/ScalarFluidInteractionFieldSlave.h	1970-01-01 00:00:00 +0000
+++ Fields/ScalarFluidInteractionFieldSlave.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,55 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+#ifndef __SCALAR_FLUID_INTERACTION_FIELD_SLAVE_H
+#define __SCALAR_FLUID_INTERACTION_FIELD_SLAVE_H
+
+// -- project includes --
+#include "Fields/FieldSlave.h"
+#include "Model/FluidInteraction.h"
+#include "pis/pi_storage_f.h"
+
+//template <typename I> class TParallelInteractionStorage;
+template <typename T> class ParallelInteractionStorage_F;
+class TML_Comm;
+
+/*!
+  \class ScalarFluidInteractionFieldSlave
+  \brief class for slave part of scalar field defined on the fluid interactions
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+
+template <typename T>
+class ScalarFluidInteractionFieldSlave : public AFieldSlave
+{
+protected:
+  virtual void SendDataWithPos();
+  virtual void SendDataWithID();
+  virtual void SendDataWithIDPos();
+  virtual void SendDataWithParticle();
+  virtual void SendDataSum();
+  virtual void SendDataMax();
+  virtual void sendData();
+
+  CFluidInteraction::ScalarFieldFunction m_rdf;
+  ParallelInteractionStorage_F<T>* m_pis;
+
+ public:
+  ScalarFluidInteractionFieldSlave(TML_Comm*,ParallelInteractionStorage_F<T>*,CFluidInteraction::ScalarFieldFunction);
+};
+
+#include "Fields/ScalarFluidInteractionFieldSlave.hpp"
+
+#endif //__SCALAR_FLUID_INTERACTION_FIELD_SLAVE_H

=== added file 'Fields/ScalarFluidInteractionFieldSlave.hpp'
--- Fields/ScalarFluidInteractionFieldSlave.hpp	1970-01-01 00:00:00 +0000
+++ Fields/ScalarFluidInteractionFieldSlave.hpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,176 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+//-- STL includes --
+#include <vector>
+#include <utility>
+
+using std::vector;
+using std::pair;
+
+// -- IO includes --
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+#include "Foundation/triplet.h"
+#include "pis/pi_storage_f.h"
+
+/*!
+  constructor
+
+  \param comm the TML communicator used for sending the data back to the master
+  \param pis a pointer to the fluid interaction storage
+  \param rdf the fluid interaction member function to access the data
+*/
+template <typename T>
+ScalarFluidInteractionFieldSlave<T>::ScalarFluidInteractionFieldSlave(TML_Comm* comm,ParallelInteractionStorage_F<T>* pis,CFluidInteraction::ScalarFieldFunction rdf):AFieldSlave(comm)
+{
+  m_rdf=rdf;
+  m_pis=pis;
+}
+
+/*!
+  Send data back to master. Determine the type of data (full/sum) to send back from
+  the received coll_type and call the according send function.
+*/
+template <typename T>
+void ScalarFluidInteractionFieldSlave<T>::sendData()
+{
+  // debug output
+  console.XDebug() << "InteractionFieldSlave<T>::sendData()\n";
+
+  int coll_type;
+  m_comm->recv_broadcast(coll_type,0);
+
+  // debug output
+  console.XDebug() << "received coll_type=" << coll_type << "\n";
+
+  switch(coll_type){
+  case COLL_TYPE_FULL_WITH_POS : SendDataWithPos();break;
+  case COLL_TYPE_FULL_WITH_ID : SendDataWithID();break;
+  case COLL_TYPE_FULL_WITH_ID_POS : SendDataWithIDPos();break;
+  case COLL_TYPE_FULL_WITH_PARTICLE : SendDataWithParticle();break;
+  case COLL_TYPE_SUM : SendDataSum();break;
+  case COLL_TYPE_MAX : SendDataMax();break;
+
+  default: std::cerr << "unknown collection type" << std::endl;
+  }
+}
+
+
+/*!
+  send full field data and position of the interaction
+*/
+template <typename T>
+void ScalarFluidInteractionFieldSlave<T>::SendDataWithPos()
+{
+  vector<pair<Vec3,double> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithPos(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+/*!
+  send full field data and id of the particle
+*/
+template <typename T>
+void ScalarFluidInteractionFieldSlave<T>::SendDataWithID()
+{
+  vector<pair<int,double> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithID(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+
+/*!
+  send full field data and id, position of the particle
+*/
+template <typename T>
+void ScalarFluidInteractionFieldSlave<T>::SendDataWithIDPos()
+{
+  vector<pair<pair<int,Vec3>,double> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithIDPos(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+/*!
+  send full field data and id, position and radius of the particle
+*/
+template <typename T>
+void ScalarFluidInteractionFieldSlave<T>::SendDataWithParticle()
+{
+  vector<pair<esys::lsm::triplet<int,Vec3,double>,double> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithParticle(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+/*!
+  send sum only
+*/
+template <typename T>
+void ScalarFluidInteractionFieldSlave<T>::SendDataSum()
+{
+  vector<double> data_vec;
+
+  // get data from interactions
+  this->m_pis->forAllInnerInteractionsGet(data_vec,m_rdf);
+
+  // sum data
+  double sum=0.0;
+  for(vector<double>::iterator iter=data_vec.begin();
+      iter!=data_vec.end();
+      iter++){
+    sum+=*iter;
+  }
+
+  vector<double> sum_vec;
+  sum_vec.push_back(sum);
+  this->m_comm->send_gather(sum_vec,0);
+}
+
+
+/*!
+  send maximum only
+*/
+template <typename T>
+void ScalarFluidInteractionFieldSlave<T>::SendDataMax()
+{
+  vector<double> data_vec;
+
+  // get data from interactions
+  this->m_pis->forAllInnerInteractionsGet(data_vec,m_rdf);
+
+  // sum data
+  double max=*(data_vec.begin());
+  for(vector<double>::iterator iter=data_vec.begin();
+      iter!=data_vec.end();
+      iter++){
+    max=(*iter > max) ? *iter : max;
+  }
+
+  vector<double> max_vec;
+  max_vec.push_back(max);
+  this->m_comm->send_gather(max_vec,0);
+}

=== added file 'Fields/VectorFluidFieldSlave.h'
--- Fields/VectorFluidFieldSlave.h	1970-01-01 00:00:00 +0000
+++ Fields/VectorFluidFieldSlave.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,47 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+#ifndef __VECTOR_FLUID_FIELD_SLAVE_H
+#define __VECTOR_FLUID_FIELD_SLAVE_H
+
+// -- project includes --
+#include "FieldSlave.h"
+#include "Model/FluidCell.h"
+
+template <class T> class ParallelParticleArray;
+class TML_Comm;
+
+/*!
+  \class VectorFluidFieldSlave
+  \brief class for slave part of vector field defined on the fluid cells
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+template <typename T>
+class VectorFluidFieldSlave : public AFieldSlave
+{
+ private:
+
+ protected: 
+  CFluidCell::VectorFieldFunction m_rdf;
+  ParallelParticleArray<T>* m_ppa;
+
+ public:
+  VectorFluidFieldSlave(TML_Comm*,ParallelParticleArray<T>*,CFluidCell::VectorFieldFunction);
+  virtual void sendData();
+};
+
+#include "VectorFluidFieldSlave.hpp"
+
+#endif //__SCALAR_FLUID_FIELD_SLAVE_H

=== added file 'Fields/VectorFluidFieldSlave.hpp'
--- Fields/VectorFluidFieldSlave.hpp	1970-01-01 00:00:00 +0000
+++ Fields/VectorFluidFieldSlave.hpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,54 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+//-- STL includes --
+#include <vector>
+#include <utility>
+
+using std::vector;
+using std::pair;
+
+// -- IO includes --
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+/*!
+  constructor 
+
+  \param comm the TML communicator used for sending the data back to the master
+  \param ppa a pointer to the particle array
+  \param rdf the fluid cell member function to access the data
+*/
+template <typename T>
+VectorFluidFieldSlave<T>::VectorFluidFieldSlave(TML_Comm* comm,ParallelParticleArray<T>* ppa,CFluidCell::VectorFieldFunction rdf):AFieldSlave(comm)
+{
+  m_ppa=ppa;
+  m_rdf=rdf;
+}
+
+/*!
+  send data back to master
+*/
+template <typename T>
+void VectorFluidFieldSlave<T>::sendData()
+{ 
+  vector<pair<Vec3,Vec3> > data_vec;
+
+  data_vec=m_ppa->forAllInnerCellsGet(m_rdf);
+ 
+  // send data to master
+  m_comm->send_gather(data_vec,0);
+}
+
+

=== added file 'Fields/VectorFluidInteractionFieldSlave.h'
--- Fields/VectorFluidInteractionFieldSlave.h	1970-01-01 00:00:00 +0000
+++ Fields/VectorFluidInteractionFieldSlave.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,54 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+#ifndef __VECTOR_FLUID_INTERACTION_FIELD_SLAVE_H
+#define __VECTOR_FLUID_INTERACTION_FIELD_SLAVE_H
+
+// -- project includes --
+#include "Fields/FieldSlave.h"
+#include "Model/FluidInteraction.h"
+#include "pis/pi_storage_f.h"
+
+template <typename T> class ParallelInteractionStorage_F;
+class TML_Comm;
+
+/*!
+  \class VectorFluidInteractionFieldSlave
+  \brief class for slave part of vector field defined on the fluid interactions
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+
+template <typename T>
+class VectorFluidInteractionFieldSlave : public AFieldSlave
+{
+protected:
+  virtual void SendDataWithPos();
+  virtual void SendDataWithID();
+  virtual void SendDataWithIDPos();
+  virtual void SendDataWithParticle();
+  virtual void SendDataSum();
+  virtual void SendDataMax();
+  virtual void sendData();
+
+  CFluidInteraction::VectorFieldFunction m_rdf;
+  ParallelInteractionStorage_F<T>* m_pis;
+
+ public:
+  VectorFluidInteractionFieldSlave(TML_Comm*,ParallelInteractionStorage_F<T>*,CFluidInteraction::VectorFieldFunction);
+};
+
+#include "Fields/VectorFluidInteractionFieldSlave.hpp"
+
+#endif //__VECTOR_FLUID_INTERACTION_FIELD_SLAVE_H

=== added file 'Fields/VectorFluidInteractionFieldSlave.hpp'
--- Fields/VectorFluidInteractionFieldSlave.hpp	1970-01-01 00:00:00 +0000
+++ Fields/VectorFluidInteractionFieldSlave.hpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,160 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+//-- STL includes --
+#include <vector>
+#include <utility>
+
+using std::vector;
+using std::pair;
+
+// -- IO includes --
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+#include "Foundation/triplet.h"
+#include "pis/pi_storage_f.h"
+
+/*!
+  constructor
+
+  \param comm the TML communicator used for sending the data back to the master
+  \param pis a pointer to the fluid interaction storage
+  \param rdf the fluid interaction member function to access the data
+*/
+template <typename T>
+VectorFluidInteractionFieldSlave<T>::VectorFluidInteractionFieldSlave(TML_Comm* comm,ParallelInteractionStorage_F<T>* pis,CFluidInteraction::VectorFieldFunction rdf):AFieldSlave(comm)
+{
+  m_rdf=rdf;
+  m_pis=pis;
+}
+
+/*!
+  Send data back to master. Determine the type of data (full/sum) to send back from
+  the received coll_type and call the according send function.
+*/
+template <typename T>
+void VectorFluidInteractionFieldSlave<T>::sendData()
+{
+  // debug output
+  console.XDebug() << "InteractionFieldSlave<T>::sendData()\n";
+
+  int coll_type;
+  m_comm->recv_broadcast(coll_type,0);
+
+  // debug output
+  console.XDebug() << "received coll_type=" << coll_type << "\n";
+
+  switch(coll_type){
+  case COLL_TYPE_FULL_WITH_POS : SendDataWithPos();break;
+  case COLL_TYPE_FULL_WITH_ID : SendDataWithID();break;
+  case COLL_TYPE_FULL_WITH_ID_POS : SendDataWithIDPos();break;
+  case COLL_TYPE_FULL_WITH_PARTICLE : SendDataWithParticle();break;
+  case COLL_TYPE_SUM : SendDataSum();break;
+  case COLL_TYPE_MAX : SendDataMax();break;
+
+  default: std::cerr << "unknown collection type" << std::endl;
+  }
+}
+
+
+/*!
+  send full field data and position of the interaction
+*/
+template <typename T>
+void VectorFluidInteractionFieldSlave<T>::SendDataWithPos()
+{
+  vector<pair<Vec3,Vec3> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithPos(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+/*!
+  send full field data and id of the particle
+*/
+template <typename T>
+void VectorFluidInteractionFieldSlave<T>::SendDataWithID()
+{
+  vector<pair<int,Vec3> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithID(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+
+/*!
+  send full field data and id, position of the particle
+*/
+template <typename T>
+void VectorFluidInteractionFieldSlave<T>::SendDataWithIDPos()
+{
+  vector<pair<pair<int,Vec3>,Vec3> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithIDPos(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+/*!
+  send full field data and id, position and radius of the particle
+*/
+template <typename T>
+void VectorFluidInteractionFieldSlave<T>::SendDataWithParticle()
+{
+  vector<pair<esys::lsm::triplet<int,Vec3,double>,Vec3> > data;
+
+  data=this->m_pis->forAllInnerInteractionsGetDataWithParticle(m_rdf);
+
+  // send data to master
+  this->m_comm->send_gather(data,0);
+}
+
+/*!
+  send sum only
+*/
+template <typename T>
+void VectorFluidInteractionFieldSlave<T>::SendDataSum()
+{
+  vector<Vec3> data_vec;
+
+  // get data from interactions
+  this->m_pis->forAllInnerInteractionsGet(data_vec,m_rdf);
+
+  // sum data
+  Vec3 sum=Vec3(0.0,0.0,0.0);
+  for(vector<Vec3>::iterator iter=data_vec.begin();
+      iter!=data_vec.end();
+      iter++){
+    sum+=*iter;
+  }
+
+  vector<Vec3> sum_vec;
+  sum_vec.push_back(sum);
+  this->m_comm->send_gather(sum_vec,0);
+}
+
+
+/*!
+  send maximum only
+*/
+template <typename T>
+void VectorFluidInteractionFieldSlave<T>::SendDataMax()
+{
+}

=== modified file 'Fields/field_const.h'
--- Fields/field_const.h	2017-01-04 08:26:43 +0000
+++ Fields/field_const.h	2019-07-28 14:55:57 +0000
@@ -23,7 +23,14 @@
     WRITE_TYPE_RAW,
     WRITE_TYPE_RAW_WITH_ID,
     WRITE_TYPE_RAW_WITH_POS_ID,
-    WRITE_TYPE_SILO
+    WRITE_TYPE_SILO,
+    /****fluid contents: begin****/
+    WRITE_TYPE_RAW_WITH_POS,
+    WRITE_TYPE_RAW_WITH_ID_POS,
+    WRITE_TYPE_RAW_WITH_PARTICLE,
+    WRITE_TYPE_VTI,
+    WRITE_TYPE_VTU
+    /****fluid contents: end****/
 };
 
 enum {
@@ -33,7 +40,12 @@
     COLL_TYPE_FULL2 = 5,
     COLL_TYPE_FULL_DX,
     COLL_TYPE_FULL_WITH_ID,
-    COLL_TYPE_FULL_WITH_POS_ID
+    COLL_TYPE_FULL_WITH_POS_ID,
+    /****fluid contents: begin****/
+    COLL_TYPE_FULL_WITH_POS,
+    COLL_TYPE_FULL_WITH_ID_POS,
+    COLL_TYPE_FULL_WITH_PARTICLE
+    /****fluid contents: end****/
 };
 
 #endif // __FIELD_CONST_H

=== added file 'Model/FluidCell.cpp'
--- Model/FluidCell.cpp	1970-01-01 00:00:00 +0000
+++ Model/FluidCell.cpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,145 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+
+// -- project includes --
+#include "Model/FluidCell.h"
+#include "tml/message/packed_message_interface.h"
+
+
+CFluidCell::CFluidCell() 
+{
+  m_Mu=0;m_Phi=0;m_newPhi=0;m_effPhi=0;m_D=0;m_K=0;m_effK=0;m_Bf=0;m_effBf=0;m_P=0;m_disP=0;m_effP=0;m_Volume=0;m_Size=Vec3(0,0,0);
+  m_Vp=Vec3(0,0,0);m_Vf=Vec3(0,0,0);m_Pg=Vec3(0,0,0);m_Pos=Vec3(0,0,0);
+  c_W=0;c_E=0;c_N=0;c_S=0;c_D=0;c_U=0;c_C=0;c_B=0;
+}
+
+
+CFluidCell& CFluidCell::operator=(const CFluidCell& rhs)
+{
+  m_Mu=rhs.m_Mu;
+  m_Phi=rhs.m_Phi;
+  m_newPhi=rhs.m_newPhi;
+  m_effPhi=rhs.m_effPhi;
+  m_D=rhs.m_D;
+  m_K=rhs.m_K;
+  m_effK=rhs.m_effK;
+  m_Bf=rhs.m_Bf;
+  m_effBf=rhs.m_effBf;
+  m_P=rhs.m_P;
+  m_disP=rhs.m_disP;
+  m_effP=rhs.m_effP;
+  m_Volume=rhs.m_Volume;
+  m_Size=rhs.m_Size;
+  m_Vp=rhs.m_Vp;
+  m_Vf=rhs.m_Vf;
+  m_Pg=rhs.m_Pg;
+  m_Pos=rhs.m_Pos;
+  m_Index=rhs.m_Index;
+  c_W=rhs.c_W;
+  c_E=rhs.c_E;
+  c_N=rhs.c_N;
+  c_S=rhs.c_S;
+  c_D=rhs.c_D;
+  c_U=rhs.c_U;
+  c_C=rhs.c_C;
+  c_B=rhs.c_B;
+  return *this;
+}
+
+
+/*!
+  Get the fluidcell member function which returns a scalar field of a given name.
+
+  \param name the name of the field
+*/
+CFluidCell::ScalarFieldFunction CFluidCell::getScalarFieldFunction(const string& name)
+{
+  CFluidCell::ScalarFieldFunction sf;
+
+  if(name=="Mu"){
+    sf=&CFluidCell::getMu;
+  } else if(name=="Phi"){
+    sf=&CFluidCell::geteffPhi;
+  } else if(name=="dPhi"){
+    sf=&CFluidCell::getdPhi;
+  } else if(name=="D"){
+    sf=&CFluidCell::getD;
+  } else if(name=="K"){
+    sf=&CFluidCell::geteffK;
+  } else if(name=="Bf"){
+    sf=&CFluidCell::geteffBf;
+  } else if(name=="P"){
+    sf=&CFluidCell::getP;
+  } else if(name=="disP"){
+    sf=&CFluidCell::getdisP;
+  } else if(name=="effP"){
+    sf=&CFluidCell::geteffP;
+  } else if(name=="dP"){
+    sf=&CFluidCell::getdP;
+  } else if(name=="Volume"){
+    sf=&CFluidCell::getVolume;
+  } else if(name=="Vp_abs"){
+    sf=&CFluidCell::getAbsVp;
+  } else if(name=="Vf_abs"){
+    sf=&CFluidCell::getAbsVf;
+  } else if(name=="c_W"){
+    sf=&CFluidCell::getc_W;
+  } else if(name=="c_E"){
+    sf=&CFluidCell::getc_E;
+  } else if(name=="c_N"){
+    sf=&CFluidCell::getc_N;
+  } else if(name=="c_S"){
+    sf=&CFluidCell::getc_S;
+  } else if(name=="c_D"){
+    sf=&CFluidCell::getc_D;
+  } else if(name=="c_U"){
+    sf=&CFluidCell::getc_U;
+  } else if(name=="c_C"){
+    sf=&CFluidCell::getc_C;
+  } else if(name=="c_B"){
+    sf=&CFluidCell::getc_B;
+  } else {
+    sf=NULL;
+    cerr << "ERROR - invalid name for fluid cell scalar  access function" << endl;
+  }
+
+  return sf;
+}
+
+
+/*!
+  Get the fluid cell member function which returns a vector field of a given name.
+
+  \param name the name of the field
+*/
+
+CFluidCell::VectorFieldFunction CFluidCell::getVectorFieldFunction(const string& name)
+  {
+  CFluidCell::VectorFieldFunction sf;
+
+  if(name=="Vp"){
+    sf=&CFluidCell::getVp;
+  } else if(name=="Vf"){
+    sf=&CFluidCell::getVf;
+  } else if(name=="Pos"){
+    sf=&CFluidCell::getPos;
+  } else if(name=="Index"){
+    sf=&CFluidCell::getIndex;
+  } else {
+    sf=NULL;
+    cerr << "ERROR - invalid name for fluid cell vector access function" << endl;
+  }
+
+  return sf;
+}
+

=== added file 'Model/FluidCell.h'
--- Model/FluidCell.h	1970-01-01 00:00:00 +0000
+++ Model/FluidCell.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,133 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+#ifndef __FLUID_CELL_H
+#define __FLUID_CELL_H
+
+// -- project includes --
+#include "Foundation/vec3.h"
+#include "tml/message/packed_message_interface.h"
+
+/*!
+  \class CFluidCell
+  \brief base class for fluid cells
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+
+class CFluidCell
+{
+ protected:
+  double m_Mu;//fluid viscosity
+  double m_Phi,m_newPhi,m_effPhi;//porosity
+  double m_D, m_K, m_effK;//permeability
+  double m_Bf, m_effBf;//fluid bulk modulus
+  double m_P, m_disP, m_effP;//pore fluid pressure
+  double m_Volume; //total particle volume
+  Vec3 m_Size;
+  Vec3 m_Vp, m_Vf;//mean particle velocity/fluid velocity
+  Vec3 m_Pg;//pressure gradient
+  Vec3 m_Pos;//global position of cell centre
+  Vec3 m_Index;//global index of cell
+
+  //coefficients for linear equations:
+  double c_W;
+  double c_E;
+  double c_N;
+  double c_S;
+  double c_D;
+  double c_U;
+  double c_C;
+  double c_B;
+
+ public:
+  CFluidCell();
+  virtual ~CFluidCell(){};
+
+  inline double getMu() const {return m_Mu;};
+  inline double getPhi() const {return m_Phi;};
+  inline double getnewPhi() const {return m_newPhi;};
+  inline double geteffPhi() const {return m_effPhi;};
+  inline double getdPhi() const {return m_newPhi-m_Phi;};
+  inline double getD() const {return m_D;};
+  inline double getK() const {return m_K;};
+  inline double geteffK() const {return m_effK;};
+  inline double getBf() const {return m_Bf;};
+  inline double geteffBf() const {return m_effBf;};
+  inline double getP() const {return m_P;};
+  inline double getdisP() const {return m_disP;};
+  inline double geteffP() const {return m_effP;};
+  inline double getdP() const {return m_disP-m_P;};
+  inline double getVolume() const {return m_Volume;};
+  inline Vec3 getSize() const {return m_Size;};
+  inline Vec3 getVp() const {return m_Vp;};
+  inline Vec3 getVf() const {return m_Vf;};
+  inline Vec3 getPg() const {return m_Pg;};
+  inline Vec3 getPos() const {return m_Pos;};
+  inline Vec3 getIndex() const {return m_Index;};
+  inline double getAbsVp() const {return m_Vp.norm();};
+  inline double getAbsVf() const {return m_Vf.norm();};
+  //get coefficients
+  inline double getc_W() const {return c_W;};
+  inline double getc_E() const {return c_E;};
+  inline double getc_N() const {return c_N;};
+  inline double getc_S() const {return c_S;};
+  inline double getc_D() const {return c_D;};
+  inline double getc_U() const {return c_U;};
+  inline double getc_C() const {return c_C;};
+  inline double getc_B() const {return c_B;};
+
+  void setMu(double Mu){m_Mu=Mu;};
+  void setPhi(double Phi){m_Phi=Phi;};
+  void setnewPhi(double newPhi){m_newPhi=newPhi;};
+  void seteffPhi(double effPhi){m_effPhi=effPhi;};
+  void setD(double D){m_D=D;};
+  void setK(double K){m_K=K;};
+  void seteffK(double effK){m_effK=effK;};
+  void setBf(double Bf){m_Bf=Bf;};
+  void seteffBf(double effBf){m_effBf=effBf;};
+  void setP(double P){m_P=P;};
+  void setdisP(double disP){m_disP=disP;};
+  void seteffP(double effP){m_effP=effP;};
+  void setVolume(double Volume){m_Volume=Volume;};
+  void setSize(Vec3 Size){m_Size=Size;};
+  void setVp(Vec3 Vp){m_Vp=Vp;};
+  void setVf(Vec3 Vf){m_Vf=Vf;};
+  void setPg(Vec3 Pg){m_Pg=Pg;};
+  void setPos(Vec3 Pos){m_Pos=Pos;};
+  void setIndex(Vec3 Index){m_Index=Index;};
+  //set coefficients
+  void setc_W(double W){c_W=W;};
+  void setc_E(double E){c_E=E;};
+  void setc_N(double N){c_N=N;};
+  void setc_S(double S){c_S=S;};
+  void setc_D(double D){c_D=D;};
+  void setc_U(double U){c_U=U;};
+  void setc_C(double C){c_C=C;};
+  void setc_B(double B){c_B=B;};
+
+  CFluidCell& operator=(const CFluidCell&);
+
+  typedef double (CFluidCell::* ScalarFieldFunction)() const;
+  typedef Vec3 (CFluidCell::* VectorFieldFunction)() const;
+
+  static ScalarFieldFunction getScalarFieldFunction(const string&);
+  static VectorFieldFunction getVectorFieldFunction(const string&);
+
+  friend class TML_PackedMessageInterface;
+
+};
+
+
+#endif //__FLUID_CELL_H

=== added file 'Model/FluidInteraction.cpp'
--- Model/FluidInteraction.cpp	1970-01-01 00:00:00 +0000
+++ Model/FluidInteraction.cpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,99 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+
+#include "Model/FluidInteraction.h"
+
+
+CFluidInteraction::CFluidInteraction (CParticle* p,CFluidCell* cell)
+{
+  m_p=p;
+  m_cell=cell;
+}
+
+
+/*!
+Calculate the fluid force.
+*/
+void CFluidInteraction::calcForces()
+{
+  double Mu=m_cell->getMu();
+  double Phi=m_cell->geteffPhi();
+  double D=m_cell->getD();
+  Vec3 vrel=m_cell->getVf()-m_p->getVel();
+  //drag force on unit volume
+  Vec3 drag;
+  if(D==0){drag=Vec3(0,0,0);}
+  else{drag=(150.0*Mu*(1.0-Phi)/Phi/D/D+1.75*1000.0*vrel.norm()/D)*vrel;};
+  //buoyancy force on unit volume
+  Vec3 buoyancy=m_cell->getPg();
+  //force applied on particle
+  double r=m_p->getRad();
+  double Volume=4.0/3.0*3.1415926*r*r*r;
+  m_drag=drag*Volume;
+  m_buoyancy=buoyancy*Volume;
+  m_force=m_drag+m_buoyancy;
+  m_p->applyForce(m_force,m_p->getPos());
+}
+
+
+/*!
+  Get the fluidinteraction member function which returns a scalar field of a given name.
+
+  \param name the name of the field
+*/
+CFluidInteraction::ScalarFieldFunction CFluidInteraction::getScalarFieldFunction(const string& name)
+{
+  CFluidInteraction::ScalarFieldFunction sf;
+  if(name=="ID"){
+    sf=&CFluidInteraction::getParticleID;
+  } else if(name=="Drag_abs"){
+    sf=&CFluidInteraction::getVbsDrag;
+  } else if(name=="Buoyancy_abs"){
+    sf=&CFluidInteraction::getVbsBuoyancy;
+  } else if(name=="Force_abs"){
+    sf=&CFluidInteraction::getVbsForce;
+  } else {
+    sf=NULL;
+    cerr << "ERROR - invalid name for fluid interaction scalar  access function" << endl;
+  }
+
+  return sf;
+}
+
+
+/*!
+  Get the fluid interaction member function which returns a vector field of a given name.
+
+  \param name the name of the field
+*/
+CFluidInteraction::VectorFieldFunction CFluidInteraction::getVectorFieldFunction(const string& name)
+{
+  CFluidInteraction::VectorFieldFunction sf;
+
+  if(name=="Position"){
+    sf=&CFluidInteraction::getParticlePos;
+  } else if(name=="Drag"){
+    sf=&CFluidInteraction::getDrag;
+  } else if(name=="Buoyancy"){
+    sf=&CFluidInteraction::getBuoyancy;
+  } else if(name=="Force"){
+    sf=&CFluidInteraction::getForce;
+  } else {
+    sf=NULL;
+    cerr << "ERROR - invalid name for fluid interaction vector access function" << endl;
+  }
+
+  return sf;
+}
+
+

=== added file 'Model/FluidInteraction.h'
--- Model/FluidInteraction.h	1970-01-01 00:00:00 +0000
+++ Model/FluidInteraction.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,72 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+
+#ifndef __FLUID_INTERACTION_H
+#define __FLUID_INTERACTION_H
+
+
+#include "Model/FluidCell.h"
+#include "Model/Particle.h"
+#include "Foundation/vec3.h"
+#include "Foundation/triplet.h"
+
+/*!
+  \class CFluidInteraction
+  \brief base class for fluid interactions
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+
+class CFluidInteraction {
+ protected:
+  CFluidCell *m_cell;
+  CParticle *m_p;
+  Vec3 m_drag;
+  Vec3 m_buoyancy;
+  Vec3 m_force;
+
+ public:
+  CFluidInteraction (CParticle*,CFluidCell*);
+  virtual ~CFluidInteraction(){};
+
+  void calcForces();
+
+  inline Vec3 getParticlePos() const {return m_p->getPos();};
+  inline Vec3 getDrag() const {return m_drag;};
+  inline Vec3 getBuoyancy() const {return m_buoyancy;};
+  inline Vec3 getForce() const {return m_force;};
+  inline double getVbsDrag() const {return m_drag.norm();};
+  inline double getVbsBuoyancy() const {return m_buoyancy.norm();};
+  inline double getVbsForce() const {return m_force.norm();};
+  inline double getParticleID() const {return m_p->getID();};
+  inline esys::lsm::triplet<int,Vec3,double> getParticleData() const
+  {
+    return
+      esys::lsm::triplet<int,Vec3,double>(
+        m_p->getID(),
+        m_p->getPos(),
+        m_p->getRad()
+      );
+  };
+
+  typedef double (CFluidInteraction::* ScalarFieldFunction)() const;
+  typedef Vec3 (CFluidInteraction::* VectorFieldFunction)() const;
+
+  static ScalarFieldFunction getScalarFieldFunction(const string&);
+  static VectorFieldFunction getVectorFieldFunction(const string&);
+
+};
+
+#endif //__FLUID_INTERACTION_H

=== modified file 'Model/Makefile.am'
--- Model/Makefile.am	2018-06-26 06:26:59 +0000
+++ Model/Makefile.am	2019-07-28 14:55:57 +0000
@@ -86,6 +86,10 @@
 	EWallInteractionGroup.hpp \
 	EWallInteraction.h \
 	EWallInteraction.hpp \
+	FluidCell.cpp \
+	FluidCell.h \
+	FluidInteraction.cpp \
+	FluidInteraction.h \
 	FractalFriction.cpp \
 	FractalFriction.h \
 	FrictionInteraction.cpp \
@@ -227,4 +231,4 @@
   $(top_builddir)/Geometry/libGgGeometry.la            \
   $(top_builddir)/tml/message/libTmlMessage.la         \
   $(top_builddir)/tml/type/libTmlType.la               \
-  $(top_builddir)/tml/comm/libTmlComm.la               
+  $(top_builddir)/tml/comm/libTmlComm.la

=== modified file 'Parallel/ASubLattice.h'
--- Parallel/ASubLattice.h	2017-09-07 08:07:19 +0000
+++ Parallel/ASubLattice.h	2019-07-28 14:55:57 +0000
@@ -46,7 +46,7 @@
   typedef std::pair<int,int>          ParticleIdPair;
   typedef std::vector<ParticleIdPair> ParticleIdPairVector;
   typedef std::vector<int> IdVector;
-  
+
   virtual ~ASubLattice();
   void setNTSize(int);
   virtual void setParticleType(const std::string &particleType)
@@ -57,6 +57,16 @@
   {
     return m_particleType;
   }
+  /****fluid contents: begin***/
+  virtual void addFluidInteraction()=0;
+  virtual void updateFluid()=0;  
+  virtual void exchangeCells()=0;
+  virtual void sendCoeffi()=0;
+  virtual void recvPressure()=0;
+  virtual void solveMatrix()=0;
+
+  /****fluid contents: end***/
+
   virtual void setTimeStepSize(double dt) = 0;
   virtual vector<int> getCommCoords() const=0;
   virtual vector<int> getCommDims() const=0;
@@ -126,8 +136,8 @@
   virtual void applyForceToWall()=0;
   virtual void setVelocityOfWall()=0;
   virtual void setParticleVelocity()=0;
-  virtual void setParticleDensity()=0; 
-  virtual void setTaggedParticleVel()=0; 
+  virtual void setParticleDensity()=0;
+  virtual void setTaggedParticleVel()=0;
   virtual void setParticleAngularVelocity(){};
   virtual void setParticleNonDynamic()=0;
   virtual void setParticleNonRot()=0;
@@ -137,7 +147,7 @@
   virtual void translateMeshBy(const std::string &meshName, const Vec3 &translation)=0;
 
   virtual void setTimer(MpiWTimers &timers) = 0;
-  
+
   virtual void do2dCalculations(bool do2d) = 0;
 
   // --- setting interaction parameters during a simulation ---
@@ -156,12 +166,18 @@
   virtual void addScalarTriangleField()=0;
   virtual void sendFieldData()=0;
   virtual void addVectorWallField()=0;
+  /****fluid contents: begin***/
+  virtual void addScalarFluidField()=0;
+  virtual void addVectorFluidField()=0;
+  virtual void addScalarFluidInteractionField()=0;
+  virtual void addVectorFluidInteractionField()=0;
+  /****fluid contents: end***/
 
   // output
   virtual void printStruct()=0;
   virtual void printData()=0;
   virtual void printTimes()=0;
-  
+
   // -- mesh data exchange --
   virtual void getMeshNodeRef()=0;
   virtual void getMeshFaceRef()=0;

=== added file 'Parallel/GMRESSolverMaster.cpp'
--- Parallel/GMRESSolverMaster.cpp	1970-01-01 00:00:00 +0000
+++ Parallel/GMRESSolverMaster.cpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,143 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+
+#include "GMRESSolverMaster.h"
+
+#include <iostream>
+#include <fstream>
+#include <math.h>
+#include <string>
+#include <sstream>
+#include <cstring>
+#include <iomanip>
+#include <algorithm>
+#include <stdlib.h>
+#include <stdio.h>
+#include "mpi.h"
+
+#include "Parallel/mpicmdbuf.h"
+#include "Parallel/mpisgvbuf.h"
+#include "Parallel/sublattice_cmd.h"
+#include "Parallel/mpibarrier.h"
+#include "Parallel/BroadCast_cmd.h"
+#include "Parallel/mpi_tag_defs.h"
+#include "tml/comm/comm.h"
+
+using namespace std;
+
+
+GMRESSolverMaster::GMRESSolverMaster(
+  MPI_Comm* comm,
+  TML_Comm* tml_comm,
+  vector<pair<Vec3,double> > co_W,
+  vector<pair<Vec3,double> > co_E,
+  vector<pair<Vec3,double> > co_N,
+  vector<pair<Vec3,double> > co_S,
+  vector<pair<Vec3,double> > co_D,
+  vector<pair<Vec3,double> > co_U,
+  vector<pair<Vec3,double> > co_C,
+  vector<pair<Vec3,double> > co_B,
+  int nx,
+  int ny,
+  int nz
+)
+{
+  m_comm=comm;
+  m_tml_comm=tml_comm;
+  MPI_Comm_rank(*m_comm, &m_rank);
+
+  Nx=nx;Ny=ny;Nz=nz;
+  Number=Nx*Ny*Nz;
+  W.resize(Number); E.resize(Number); N.resize(Number); S.resize(Number); D.resize(Number); U.resize(Number); C.resize(Number); B.resize(Number);
+  irow.resize(Number*7); icol.resize(Number*7); var.resize(Number*7);
+  index.resize(Number);
+  for(int n=0;n<Number;n++){
+    W[n]=(&co_W[n])->second; E[n]=(&co_E[n])->second;
+    N[n]=(&co_N[n])->second; S[n]=(&co_S[n])->second;
+    D[n]=(&co_D[n])->second; U[n]=(&co_U[n])->second;
+    C[n]=(&co_C[n])->second; B[n]=(&co_B[n])->second;
+    index[n]=(&co_W[n])->first;
+  }
+}
+
+
+GMRESSolverMaster::~GMRESSolverMaster()
+{
+  irow.clear();icol.clear();
+  W.clear();E.clear();N.clear();S.clear();D.clear();U.clear();C.clear();B.clear();
+  var.clear();index.clear();
+}
+
+
+vector<pair<Vec3,double> > GMRESSolverMaster::MatrixSolving()
+{
+  CMPILCmdBuffer cmd_buffer(*m_comm,m_rank);
+  CMPIBarrier barrier(*m_comm);
+
+  GlobalMatrixSetting();
+
+  cmd_buffer.broadcast(CMD_SOLVE);
+
+  m_tml_comm->broadcast_cont(B);
+  m_tml_comm->broadcast_cont(irow);
+  m_tml_comm->broadcast_cont(icol);
+  m_tml_comm->broadcast_cont(var);
+  m_tml_comm->broadcast(Number);
+  m_tml_comm->broadcast(nnz);
+
+  multimap<int,double> temp;
+  m_tml_comm->gather(temp);
+
+  barrier.wait("solve matrix");
+
+  vector<double> x;
+  for(multimap<int,double>::iterator iter=temp.begin();
+      iter!=temp.end();
+      iter++){
+    x.push_back(iter->second);
+  };
+
+  //putting x into results vector
+  vector<pair<Vec3,double> > results;
+  for(int n=0;n<Number;n++){
+    results.push_back(make_pair(index[n],x[n]));
+  }
+
+  x.clear();
+  temp.clear();
+  return results;
+}
+
+
+void GMRESSolverMaster::GlobalMatrixSetting()
+{
+  nnz=0;
+  for(int n=1,m1=0,m2=0;n<=Number;n++){
+    if (n>Nx*Ny){irow[nnz]=n;icol[nnz]=n-Nx*Ny;var[nnz]=D[n-1];nnz++;};//down3
+    if (n>Nx*(1+Ny*m1) && n<=Nx*Ny*(m1+1)) {
+      irow[nnz]=n;icol[nnz]=n-Nx;var[nnz]=N[n-1];nnz++; //down2
+      if(n==Nx*Ny*(m1+1)){m1++;};
+    };
+    if(n%Nx!=1){irow[nnz]=n;icol[nnz]=n-1;var[nnz]=W[n-1];nnz++;};//down1
+    irow[nnz]=n;icol[nnz]=n;var[nnz]=C[n-1];nnz++; //centre
+    if(n%Nx!=0){irow[nnz]=n;icol[nnz]=n+1;var[nnz]=E[n-1];nnz++;};//up1
+    if (n>m2*Nx*Ny && n<=(m2+1)*Ny*Nx-Nx) {
+      irow[nnz]=n;icol[nnz]=n+Nx;var[nnz]=S[n-1];nnz++;//up2
+      if(n==(m2+1)*Ny*Nx-Nx){m2++;};
+    };
+    if (n<=Nx*Ny*(Nz-1)){irow[nnz]=n;icol[nnz]=n+Nx*Ny;var[nnz]=U[n-1];nnz++;};//up3
+  };
+}
+
+
+

=== added file 'Parallel/GMRESSolverMaster.h'
--- Parallel/GMRESSolverMaster.h	1970-01-01 00:00:00 +0000
+++ Parallel/GMRESSolverMaster.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,81 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+
+#ifndef __GMRESSOLVERMASTER_H
+#define __GMRESSOLVERMASTER_H
+
+//--- STL includes ---
+#include <string>
+using std::string;
+
+// --- project includes ---
+#include "Foundation/console.h"
+#include "Foundation/vec3.h"
+
+// --- TML includes ---
+#include "tml/comm/comm.h"
+
+class TML_Comm;
+
+/*!
+  \class GMRESSolverMaster
+  \brief Master class of GMRES Solver for solving sparse matrix of linear equations
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+class GMRESSolverMaster
+{
+ private:
+  int Nx,Ny,Nz,Number,nnz;
+  vector<int> irow,icol;
+  vector<double> W, E, N, S, D, U, C, B;
+  vector <double> var;
+  vector<Vec3> index;
+  int m_rank;
+
+ protected:
+  MPI_Comm *m_comm;
+  TML_Comm *m_tml_comm;
+
+  void GlobalMatrixSetting();
+
+ public:
+  vector<pair<Vec3,double> > MatrixSolving();
+
+  GMRESSolverMaster(
+    MPI_Comm*,
+    TML_Comm*,
+    vector<pair<Vec3,double> >,
+    vector<pair<Vec3,double> >,
+    vector<pair<Vec3,double> >,
+    vector<pair<Vec3,double> >,
+    vector<pair<Vec3,double> >,
+    vector<pair<Vec3,double> >,
+    vector<pair<Vec3,double> >,
+    vector<pair<Vec3,double> >,
+    int,
+    int,
+    int
+  );
+  ~GMRESSolverMaster();
+};
+
+#endif  //__GMRESSOLVERMASTER_H
+
+
+
+
+
+

=== added file 'Parallel/GMRESSolverSlave.cpp'
--- Parallel/GMRESSolverSlave.cpp	1970-01-01 00:00:00 +0000
+++ Parallel/GMRESSolverSlave.cpp	2019-07-28 14:55:57 +0000
@@ -0,0 +1,573 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+
+
+#include "GMRESSolverSlave.h"
+
+#include <iostream>
+#include <fstream>
+#include <math.h>
+#include <string>
+#include <sstream>
+#include <cstring>
+#include <iomanip>
+#include <algorithm>
+#include <stdlib.h>
+#include <stdio.h>
+#include "mpi.h"
+
+using namespace std;
+
+
+
+GMRESSolverSlave::GMRESSolverSlave(
+  MPI_Comm* comm,
+  TML_Comm* tml_comm,
+  vector<double> B0,
+  vector<int> row,
+  vector<int> col,
+  vector<double> v,
+  int number,
+  int nz
+)
+{
+  m_worker_comm=comm;
+  m_tml_global_comm=tml_comm;
+  MPI_Comm_rank(*m_worker_comm, &id);
+  MPI_Comm_size(*m_worker_comm, &nproc);
+
+  Number=number;
+  nrowsI=Number/nproc;
+  B.resize(Number);
+  for(int n=0;n<Number;n++){
+    B[n]=B0[n];
+  }
+  nnz=nz;
+  irow.resize(nnz); icol.resize(nnz);var.resize(nnz);
+  for(int n=0;n<nnz;n++){
+    irow[n]=row[n];
+    icol[n]=col[n];
+    var[n]=v[n];
+  }
+}
+
+GMRESSolverSlave::~GMRESSolverSlave()
+{
+  vector<int> zero;
+  irow.swap(zero);icol.swap(zero);ia.swap(zero);ja.swap(zero);
+  vector<double> zero_db;
+  var.swap(zero_db);A.swap(zero_db);B.swap(zero_db);b.swap(zero_db);
+}
+
+
+
+void GMRESSolverSlave::LocalMatrixSolving()
+{
+  /* local vars */
+  int precond_flag = 0;//0 for normal GMRES and 1 for preconditioned GMRES algorithm
+  double gmres_res = 1.e-5; //gmres residual tolerance
+  int i = 0, j = 0, k, ss,  k_end = 200, mm = 10, k_sum=0;//k_end=number of GMRES iterations; m=number of restarts;
+
+  LocalMatrixSetting();
+
+  //localize
+  int **receive = NULL;
+  int *cntowner= NULL;
+  int **to_be_sent= NULL;
+  int *cntsent = NULL;
+  localize(&receive, &cntowner, &to_be_sent, &cntsent);
+
+  int n_total = max_int_array(ja, my_nnz)+1;
+  double norm_r0 = 0.;
+
+  double *u = (double *)calloc( nrowsI , sizeof(double) );
+  double **v = (double **)malloc( 1 * sizeof(double *) );
+  double **h = (double **)malloc( 1 * sizeof(double *) );
+  double *x0 = (double *)calloc( n_total , sizeof(double) );
+  double *r0 = (double *)calloc( nrowsI , sizeof(double) );
+  double *g = (double *)calloc( 1 , sizeof(double) );
+  double *c = NULL, *s = NULL, *resi = NULL;
+  double delta = 0., gamma = 0.;
+
+  vector<double> x_local;
+  x_local.resize(nrowsI);
+
+  //pick the diagonal entries for preconditioning (if it is preconditioned)
+  double *diag = (double *) calloc( nrowsI , sizeof(double) );
+  double *yk = (double *) calloc( n_total , sizeof(double) );
+  if( precond_flag ) //if it is preconditioned
+  {
+    for (i = 0 ; i < nrowsI; i++ )
+      for (j = ia[i] ; j < ia[i+1]; j++)
+        if(ja[j] == i )
+          diag[i] = A[j];
+  }
+
+  //init x0
+  for(i = 0; i < nrowsI; i++)
+    x0[i] = 1.0;
+
+  //restarting vars
+  int rst = 0,  completed = 0, total_itr = 1;
+
+  //restarting loop
+  for(rst = 0; rst < mm; rst++){
+    gather_sync(receive, cntowner, to_be_sent, cntsent, x0);
+    //matrix vector product
+    sparse_matmult(x0, r0);
+    //updating the final r0
+    for(i = 0; i < nrowsI; i++)
+      r0[i] = b[i] - r0[i];
+
+    // normalization ...
+    dotproduct(r0, r0, &norm_r0); //dot product
+    norm_r0 = sqrt(norm_r0); //root
+    for(i = 0; i < nrowsI; i++)
+      r0[i] /= norm_r0; //normalizing
+
+    // initial allocation of v
+    h[0] = (double *)calloc(1 , sizeof(double) );
+    v[0] = (double *)calloc(n_total , sizeof(double) );
+    for(i = 0; i < nrowsI; i++)
+      v[0][i] = r0[i]; //initializing ...
+
+    g[0] = norm_r0; //initial g
+
+    //gmres loop
+    for(k = 0; k < k_end; k++)
+    {
+      //reallocating vars to have enough space for the next items
+      g = (double *) realloc (g, (k+2) * sizeof(double) );
+      g[k+1] = 0.;
+      c = (double *) realloc (c, (k+1) * sizeof(double) );
+      s = (double *) realloc (s, (k+1) * sizeof(double) );
+      resi = (double *) realloc (resi, (k+1) * sizeof(double) );
+      v = (double **)realloc (v, (k+2) * sizeof(double *) );
+      v[k+1] = (double *)calloc(n_total , sizeof(double) );
+      h = (double **)realloc (h , (k+2) * sizeof(double *) );
+      h[k+1] = (double *)calloc(1 , sizeof(double) );
+
+      for (j = 0; j <= (k+1); j++)
+        h[j] = (double *)realloc( h[j], (k+1) *sizeof(double) );
+
+      if( precond_flag )
+      {
+        for (i = 0 ; i < nrowsI; i++ )
+          if( diag[i] )
+            yk[i] = v[k][i] / diag[i];
+          else
+            yk[i] = v[k][i];
+        //gather/sync
+        gather_sync(receive, cntowner, to_be_sent, cntsent, yk);
+        //matrix vector product
+        sparse_matmult(yk, u);
+      }
+      else
+      {
+        //gather/sync
+        gather_sync(receive, cntowner, to_be_sent, cntsent, v[k]);
+        //matrix vector product
+        sparse_matmult(v[k], u);
+      }
+
+      for (j = 0 ; j <= k; j++)
+      {
+        dotproduct(v[j], u, &h[j][k]); //dot product
+        for( ss = 0; ss < nrowsI; ss++)
+          u[ss] -= h[j][k] * v[j][ss];
+      }
+      dotproduct(u, u, &h[k+1][k]); //dot product
+      h[k+1][k] = sqrt(h[k+1][k]); //norm
+      //updating v[k+1]
+      for(ss = 0; ss < nrowsI; ss++)
+        v[k+1][ss] = u[ss] / h[k+1][k];
+
+      for (j = 0 ; j < k; j++)
+      {
+        delta = h[j][k];
+        h[j][k] = c[j] * delta + s[j] * h[j+1][k];
+        h[j+1][k] = -s[j] * delta + c[j] * h[j+1][k];
+      }
+      gamma = sqrt(h[k][k] * h[k][k] + h[k+1][k]*h[k+1][k]);
+      c[k] = h[k][k] / gamma;
+      s[k] = h[k+1][k] / gamma;
+      h[k][k] = gamma;
+      h[k+1][k] = 0.;
+      delta = g[k];
+      g[k] = c[k] * delta + s[k] * g[k+1];
+      g[k+1] = -s[k] * delta + c[k] * g[k+1];
+      resi[k] = fabs(g[k+1]);
+
+      if (resi[k] <= gmres_res)
+      {
+        completed = 1; //set the completed flag
+        break; //terminate gmres loop ASAP!
+      }
+
+      if (rst==mm-1 && k==k_end-1 && resi[k] > gmres_res && id==0){
+        cerr << "Fails to converge, please try a larger number of GMRES iterations or restarts \n";
+        exit(0);
+      };
+
+    } //end of the main gmres loop
+
+    k--; //reducing k to emulate the inside the loop effect
+    k_sum+=k;
+
+    //compute alpha
+    double *alpha = (double *)calloc(k+1 , sizeof(double) );
+    //solve backward
+    for (j = k ; j >= 0; j--)
+    {
+      alpha[j] = g[j]/h[j][j];
+      for (ss = (j+1) ; ss <= k; ss++)
+        alpha[j] -= (h[j][ss]/h[j][j] * alpha[ss]);
+    }
+    //compute zk
+    double *zk = (double *)calloc(nrowsI , sizeof(double) );
+    for (j = 0 ; j <= k; j++)
+      for (ss = 0 ; ss < nrowsI; ss++)
+        zk[ss] += alpha[j] * v[j][ss];
+
+    if(precond_flag )
+      for (i = 0 ; i < nrowsI; i++ )
+        if (diag[i] )
+          zk[i] /= diag[i];
+
+    //compute solution
+    for (ss = 0 ; ss < nrowsI; ss++){
+      x_local[ss] = x0[ss] + zk[ss];
+    }
+
+    //preparing for restart operation ...
+    for(i = 0; i < nrowsI; i++)
+	  x0[i] = x_local[i]; //put last x_local into x0
+
+    //clean ups
+    free(alpha);
+    free(zk);
+    for (j = 0 ; j <= (k+1); j++)
+    {
+      free(v[j]);
+      free(h[j]);
+    }
+
+    if(completed)
+      break; //terminate restart loop ASAP!
+
+  } //end of restarting loop
+  //Sending local results for gather by master
+  m_tml_global_comm->send_gather(x_local,0);
+
+  free(u);
+  free(x0);
+  free(r0);
+  free(g);
+  free(c);
+  free(diag);
+  free(yk);
+  free(s);
+  free(resi);
+  free(v);
+  free(h);
+  free(cntowner);
+  free(cntsent);
+  for (j = 0 ; j < nproc; j++)
+  {
+      free(receive[j]);
+      free(to_be_sent[j]);
+  }
+  free(receive);
+  free(to_be_sent);
+}
+
+
+
+void GMRESSolverSlave::LocalMatrixSetting()
+{
+  //Setting local CRS matrix
+  int i, i_start, i_end, DI;
+
+  DI = Number/nproc;
+  int *col_at_row = (int *)calloc( DI , sizeof(int) );
+  i_start = id * DI;
+  i_end = (id+1) * DI;
+
+  int irow2,icol2;
+  for(i = 0; i < nnz; i++) {
+    irow2=irow[i]-1; //converting to zero-based
+    if ( (i_start <= irow2) && (irow2 < i_end) ) {
+      col_at_row[(irow2-i_start)]++;
+    };
+  }
+
+  //find my total number of non-zero items (my_nnz) and also update ia[]
+  ia.resize(DI+1);
+  my_nnz=0;
+  for(i= 0; i < DI; i++){
+    ia[i] = my_nnz;
+    my_nnz += col_at_row[i];
+  }
+  ia[i] = my_nnz;
+
+  //allocate A one big chunck of doubles for local CRS array for storing data
+  A.resize(my_nnz);
+  //allocate (*ja)
+  ja.resize(my_nnz);
+
+  //fill data matrix A[], ia[] and ja[]
+  for(i= 0; i < DI; i++)
+    col_at_row[i] = 0; //reset col_at_row[]
+
+  for( i = 0; i < nnz; i++){
+    irow2=irow[i]-1;
+    icol2=icol[i]-1;//converting to zero-based
+    if ((i_start <= irow2) && (irow2 < i_end) ){
+      irow2 -= i_start;
+      A[ia[irow2] + col_at_row[irow2]] = var[i];
+      ja[ia[irow2] + col_at_row[irow2]] = icol2;
+      col_at_row[irow2]++;
+    }
+  }
+
+  //setting local rhs
+  b.resize(nrowsI);
+  for(i = 0; i < Number; i++) {
+    if ( (i_start <= i) && (i < i_end) )
+      b[(i-i_start)] = B[i]; //put it in local [b]
+  }
+
+  //clean-up
+  free(col_at_row);
+}
+
+
+//finds the maximum of the given integer array with size "n"
+int GMRESSolverSlave::max_int_array(vector<int> input, int n)
+{
+  //locals
+  int i = 0;
+  int max = input[0];
+  for( i = 0; i < n ; i++)
+    if( input[i] > max)
+      max = input[i];
+  //returning the maximum value
+  return max;
+}
+
+
+void GMRESSolverSlave::localize(int ***receive, int **cntowner, int ***to_be_sent, int **cntsent)
+{
+  // locals
+  int i, j, indx, owner;
+  int phn = nrowsI;
+  int *tag = (int *)calloc(Number , sizeof(int) );
+  (*receive) = (int **)malloc(nproc * sizeof(int *));
+  (*to_be_sent) = (int **)malloc(nproc * sizeof(int *));
+  int **expect = (int **)malloc(nproc * sizeof(int *));
+  (*cntowner) = (int *)calloc(nproc , sizeof(int));
+  (*cntsent) = (int *)calloc(nproc , sizeof(int));
+
+  int *tmp_receive = NULL, *tmp_expect = NULL;
+  int msg_tag1 = 40, msg_tag2 = 60;
+  MPI_Status status1[nproc];
+  MPI_Request request1[nproc];
+  MPI_Status status2[nproc];
+  MPI_Request request2[nproc];
+
+  //initializing to NULL (safe)
+  for ( i = 0; i < nproc; i++){
+    (*receive)[i] = NULL;
+    expect[i] = NULL;
+    (*to_be_sent)[i] = NULL;
+  }
+
+  // main loop for hacking ja[]
+  for ( i = 0; i < nrowsI; i++) { // i-local always starts from zero
+    for ( indx = ia[i]; indx < ia[i+1]; indx++) //ia[] always starts from zero for this process
+    {
+      j=ja[indx]; //find the global column (zero-based)
+      owner = (int)(j/nrowsI); //and the owner rank
+      if (owner == id)
+        ja[indx] -= (id*nrowsI); // NOTE: my_rank*nlocal=vect_start_indx for this process
+      else
+      {
+        if( !tag[j] ) //is it already localized? if no do it!
+        {
+	  //updating receive list
+  	  tmp_receive = (int *) realloc ((*receive)[owner],  ((*cntowner)[owner]+1)*sizeof(int));
+   	  if(tmp_receive == NULL) {
+	    cout<<"can't realloc receive list on processor "<<id<<endl;
+	    fflush(stdout);
+	    exit(0);
+	  } else
+	    (*receive)[owner] = tmp_receive;
+
+	  (*receive)[owner][(*cntowner)[owner]] = phn;
+	  //updating expect list
+	  tmp_expect = (int *) realloc (expect[owner],  ((*cntowner)[owner]+1) * sizeof(int));
+	  if(tmp_expect == NULL) {
+            cout<<"can't realloc expect list on processor "<<id<<endl;
+	    fflush(stdout);
+	    exit(0);
+	  } else
+	    expect[owner] = tmp_expect;
+
+	  expect[owner][(*cntowner)[owner]] = j;
+	  // updating countors
+	  tag[j] = phn;
+	  ((*cntowner)[owner])++;
+	  phn++;
+	}
+        ja[indx] = tag[j];
+      }
+    }
+  }
+
+  // sending expect list to other processes
+  int rq_inc1 = 0, rq_inc2 = 0; //increaments for requests
+  for( owner = 0; owner < nproc; owner++) //loop over process which we want to send expect list
+    if ( owner == id) //dont send to myself
+      continue;
+    else
+    {
+      //sending the numbers of expected values to that process
+      MPI_Isend(&(*cntowner)[owner], 1, MPI_INT, owner, msg_tag1, *m_worker_comm, request1+rq_inc1++);
+      //sending expected values to that process
+      MPI_Isend(expect[owner], (*cntowner)[owner], MPI_INT, owner, msg_tag2, *m_worker_comm, request2+ rq_inc2++);
+    }
+
+  // receiving to_be_sent list from other processes
+  int st_inc1 = 0, st_inc2 = 0; //increaments for status
+  for( owner = 0; owner < nproc; owner++) //loop over process which we want to receive from them.
+    if ( owner == id)
+      continue;
+    else
+    {
+      //receiving the numbers of expected values for that process
+      MPI_Recv(&(*cntsent)[owner], 1, MPI_INT, owner, msg_tag1, *m_worker_comm, status1+st_inc1++ );
+      //once get the size of data, reallocating enough space
+      tmp_expect = (int *)realloc ((*to_be_sent)[owner], (*cntsent)[owner]*sizeof(int) );
+      if(tmp_expect == NULL){
+        printf("\ncan't realloc to_be_sent list!\n");
+        exit(0);
+      }
+      else {
+        (*to_be_sent)[owner] = tmp_expect;
+      }
+      // reciving data and putting them in place
+      MPI_Recv((*to_be_sent)[owner],(*cntsent)[owner], MPI_INT, owner, msg_tag2, *m_worker_comm, status2+st_inc2++ );
+      for ( i = 0; i < (*cntsent)[owner]; i++) //localizing data
+        (*to_be_sent)[owner][i] -= (id*nrowsI);
+    }
+
+  // wait until all send and receives ar complete
+  MPI_Waitall(nproc-1, request1, status1);
+  MPI_Waitall(nproc-1, request2, status2);
+
+  // clean-up
+  for( i = 0; i < nproc; i++)
+    if( i != id)
+      free(expect[i]);
+  free(expect);
+  free(tag);
+}
+
+
+
+// for gathering data in one location in the local process and synchronization among processes
+void GMRESSolverSlave::gather_sync(int **receive, int *cntowner, int **to_be_sent, int *cntsent, double *qk)
+{
+  int i = 0;
+  int owner = 0;
+  int msg_tag = 40;
+  MPI_Status status[nproc];
+  MPI_Request request[nproc];
+  double **buff_snd = (double **)malloc(nproc * sizeof(double *));
+  double **buff_rcv = (double **)malloc(nproc * sizeof(double *));
+
+  // sending doubles to other processes
+  int rq_inc = 0; //increament for requests
+  for(owner = 0; owner < nproc; owner++) //loop over process which we want to send
+    if (owner == id) //dont send to myself
+      continue;
+    else
+    {
+      // first pack doubles into the send buffer
+      buff_snd[owner] = (double *)malloc(cntsent[owner] * sizeof(double));
+      for( i = 0; i < cntsent[owner]; i++)
+	    buff_snd[owner][i] = qk[to_be_sent[owner][i]];
+
+      //sending expected values to that process
+      MPI_Isend(buff_snd[owner], cntsent[owner], MPI_DOUBLE, owner, msg_tag, *m_worker_comm, request+rq_inc++);
+    }
+
+  // receiving doubles from other processes and putting them in the right places
+  int st_inc = 0; //increament for status
+  for( owner = 0; owner < nproc; owner++) //loop over process which we want to receive from them.
+    if ( owner == id)
+      continue;
+    else
+    {
+      // first allocate receive buffer
+      buff_rcv[owner] = (double *)malloc(cntowner[owner] * sizeof(double));
+      // reciving data and putting them in receive buffer
+      MPI_Recv(buff_rcv[owner], cntowner[owner], MPI_DOUBLE, owner, msg_tag, *m_worker_comm, status+st_inc++ );
+      for ( i = 0; i < cntowner[owner]; i++) //putting in right place indicated by ja[]
+        qk[receive[owner][i]] = buff_rcv[owner][i];
+    }
+
+  // wait until all send and receives ar complete
+  MPI_Waitall(nproc-1, request, status);
+
+  //clean - ups
+  for( i = 0; i < nproc ; i++)
+    if( i != id) {
+      free(buff_snd[i]);
+      free(buff_rcv[i]);
+    }
+    free(buff_rcv);
+    free(buff_snd);
+}
+
+
+void GMRESSolverSlave::sparse_matmult(double *q, double *u)
+{
+  //locals
+  int r, j;
+  // reset u
+  for( r = 0 ; r < nrowsI; r++)
+    u[r] = 0.;
+
+  for( r = 0 ; r < nrowsI; r++)
+    for ( j = ia[r]; j < ia[r+1]; j++) {
+     u[r] += (A[j] * q[ja[j]]);
+    }
+}
+
+
+void GMRESSolverSlave::dotproduct(double *a, double *b, double *c)
+{
+  //locals
+  int i = 0;
+  double local_dot = 0.;
+  //perform local dot product
+  for( i = 0; i < nrowsI; i++)
+    local_dot += (a[i] * b[i]);
+
+  //reduceall to single values
+  MPI_Allreduce(&local_dot,c,1,MPI_DOUBLE,MPI_SUM,*m_worker_comm);
+}
+
+
+

=== added file 'Parallel/GMRESSolverSlave.h'
--- Parallel/GMRESSolverSlave.h	1970-01-01 00:00:00 +0000
+++ Parallel/GMRESSolverSlave.h	2019-07-28 14:55:57 +0000
@@ -0,0 +1,76 @@
+/////////////////////////////////////////////////////////////
+//                                                         //
+// Copyright (c) 2003-2014 by The University of Queensland //
+// Centre for Geoscience Computing                         //
+// http://earth.uq.edu.au/centre-geoscience-computing      //
+//                                                         //
+// Primary Business: Brisbane, Queensland, Australia       //
+// Licensed under the Open Software License version 3.0    //
+// http://www.apache.org/licenses/LICENSE-2.0          //
+//                                                         //
+/////////////////////////////////////////////////////////////
+
+
+#ifndef __GMRESSOLVERSLAVE_H
+#define __GMRESSOLVERSLAVE_H
+
+//--- STL includes ---
+#include <string>
+using std::string;
+
+// --- project includes ---
+#include "Foundation/console.h"
+#include "Foundation/vec3.h"
+
+// --- TML includes ---
+#include "tml/comm/comm.h"
+
+class TML_Comm;
+
+/*!
+  \class GMRESSolverSlave
+  \brief Slave class of GMRES Solver for solving sparse matrix of linear equations
+
+  \author Qi Shao
+  $Revision$
+  $Date$
+*/
+class GMRESSolverSlave
+{
+ private:
+  int nrowsI,nproc,id,nnz;
+  int Nx,Ny,Nz,Number,my_nnz;
+  vector<int> irow,icol,ia,ja;
+  vector<double> var, A, B, b;
+
+ protected:
+  MPI_Comm *m_worker_comm;
+  TML_Comm *m_tml_global_comm;
+
+  void LocalMatrixSetting();
+  int max_int_array(vector<int>, int);
+  void localize(int ***, int **, int ***, int **);
+  void gather_sync(int **, int *, int **, int *, double *);
+  void sparse_matmult(double *, double *);
+  void dotproduct(double *, double *, double *);
+
+ public:
+  void LocalMatrixSolving();
+
+  GMRESSolverSlave(
+    MPI_Comm*,
+    TML_Comm*,
+    vector<double>,
+    vector<int>,
+    vector<int>,
+    vector<double>,
+    int,
+    int
+  );
+   ~GMRESSolverSlave();
+};
+
+#endif //__GMRESSOLVERSLAVE_H
+
+
+

=== modified file 'Parallel/LatticeMaster.cpp'
--- Parallel/LatticeMaster.cpp	2019-03-12 16:42:13 +0000
+++ Parallel/LatticeMaster.cpp	2019-07-28 14:55:57 +0000
@@ -97,8 +97,12 @@
     m_center_id(0),
     m_total_time(0),
     m_t(0),
+    m_t_f(0), //fluid contents
+    m_t_old(0), //fluid contents
     m_dt(0),
+    m_dt_f(0), //fluid contents
     m_isInitialized(false),
+    m_fluidinitiated(false), //fluid contents
     m_particle_type(),
     m_preRunnableVector(),
     m_postRunnableVector(),
@@ -133,6 +137,291 @@
   console.Debug() << "CLatticeMaster::~CLatticeMaster: exit\n";
 }
 
+/****fluid contents: begin****/
+/*!
+  add fluid to the sublattice
+
+  \param Bw Bulk modulus of water
+  \param Bp Bulk modulus of particle
+  \param Mu Viscosity of water
+  \param alpha Adjusting factor between two time steps
+*/
+void CLatticeMaster::addFluidInteraction(double cellside,double Bw, double Bp, double Mu, double alpha, double flowrate, double pressure, Vec3 inflow, Vec3 outflow, double dt_f)
+{
+  console.Debug() << "CLatticeMaster::addFluidInteraction\n";
+
+  BroadcastCommand cmd(getGlobalRankAndComm(), CMD_ADDFLUID);
+
+  cmd.append(cellside);
+  cmd.append(Bw);
+  cmd.append(Bp);
+  cmd.append(Mu);
+  cmd.append(alpha);
+  cmd.append(flowrate);
+  cmd.append(pressure);
+  cmd.append(inflow);
+  cmd.append(outflow);
+  //send command to slave
+  cmd.broadcast();
+  m_dt_f=dt_f;
+  m_fluidinitiated=true;
+  console.Debug() << "end CLatticeMaster::addFluidInteraction() \n";
+}
+
+/*!
+  add fluid to the sublattice
+
+  \param Bw Bulk modulus of water
+  \param Bp Bulk modulus of particle
+  \param Mu Viscosity of water
+  \param alpha Adjusting factor between two time steps
+*/
+void CLatticeMaster::addFluidInteractionVec3(Vec3 cellside,double Bw, double Bp, double Mu, double alpha, double flowrate, double pressure, Vec3 inflow, Vec3 outflow, double dt_f)
+{
+  console.Debug() << "CLatticeMaster::addFluidInteractionVec3\n";
+
+  BroadcastCommand cmd(getGlobalRankAndComm(), CMD_ADDFLUID_VEC3);
+
+  cmd.append(cellside);
+  cmd.append(Bw);
+  cmd.append(Bp);
+  cmd.append(Mu);
+  cmd.append(alpha);
+  cmd.append(flowrate);
+  cmd.append(pressure);
+  cmd.append(inflow);
+  cmd.append(outflow);
+  //send command to slave
+  cmd.broadcast();
+  m_dt_f=dt_f;
+  m_fluidinitiated=true;
+  console.Debug() << "end CLatticeMaster::addFluidInteractionVec3() \n";
+}
+
+
+/*!
+  add a scalar fluid field to the list of fields to be saved
+
+  \param filename the name of the file the field is saved into
+  \param fieldname the name of the field
+  \param savetype output file format
+  \param t_0 first timestep to be saved
+  \param t_end last timestep to be saved
+  \param dt timesteps between saves
+*/
+void CLatticeMaster::addScalarFluidSaveField(
+  const string& filename,
+  const string& fieldname,
+  const string& savetype,
+  int t_0,
+  int t_end,
+  int dt
+)
+{
+  console.Debug()
+    << "CLatticeMaster::addScalarFluidSaveField("
+    << filename
+    << ","
+    << fieldname
+    << ","
+    << t_0
+    << ","
+    << t_end
+    << ","
+    << dt
+    << ")\n";
+
+  CMPILCmdBuffer cmd_buffer(m_global_comm,m_global_rank);
+  CMPIBarrier barrier(m_global_comm);
+
+  //send cmd to slave
+  cmd_buffer.broadcast(CMD_ADD_SFF);
+
+  AFieldMaster* new_fm =
+    new ScalarFluidFieldMaster(
+      &m_tml_global_comm,
+      fieldname,
+      filename,
+      savetype,
+      t_0,
+      t_end,
+      dt
+    );
+  m_save_fields.push_back(new_fm);
+
+  barrier.wait("CLatticeMaster::addScalarFluidSaveField");
+
+  console.Debug() << "end CLatticeMaster::addScalarFluidSaveField()\n";
+}
+
+/*!
+  add a vector fluid field to the list of fields to be saved
+
+  \param filename the name of the file the field is saved into
+  \param fieldname the name of the field
+  \param savetype output file format
+  \param t_0 first timestep to be saved
+  \param t_end last timestep to be saved
+  \param dt timesteps between saves
+*/
+void CLatticeMaster::addVectorFluidSaveField(
+  const string& filename,
+  const string& fieldname,
+  const string& savetype,
+  int t_0,
+  int t_end,
+  int dt
+)
+{
+  console.Debug()
+    << "CLatticeMaster::addVectorFluidSaveField("
+    << filename
+    << ","
+    << fieldname
+    << ","
+    << t_0
+    << ","
+    << t_end
+    << ","
+    << dt
+    << ")\n";
+
+  CMPILCmdBuffer cmd_buffer(m_global_comm,m_global_rank);
+  CMPIBarrier barrier(m_global_comm);
+
+  //send cmd to slave
+  cmd_buffer.broadcast(CMD_ADD_VFF);
+
+  AFieldMaster* new_fm =
+    new VectorFluidFieldMaster(
+      &m_tml_global_comm,
+      fieldname,
+      filename,
+      savetype,
+      t_0,
+      t_end,
+      dt
+    );
+  m_save_fields.push_back(new_fm);
+
+  barrier.wait("CLatticeMaster::addVectorFluidSaveField");
+
+  console.Debug() << "end CLatticeMaster::addVectorFluidSaveField()\n";
+}
+
+/*!
+  add a scalar fluid interaction field to the list of fields to be saved
+
+  \param filename the name of the file the field is saved into
+  \param fieldname the name of the field
+  \param savetype output file format
+  \param t_0 first timestep to be saved
+  \param t_end last timestep to be saved
+  \param dt timesteps between saves
+*/
+void CLatticeMaster::addScalarFluidInteractionSaveField(
+  const string& filename,
+  const string& fieldname,
+  const string& savetype,
+  int t_0,
+  int t_end,
+  int dt
+)
+{
+  console.Debug()
+     << "CLatticeMaster::addScalarFluidInteractionSaveField("
+     << filename
+     << ","
+     << fieldname
+     << ","
+     << t_0
+     << ","
+     << t_end
+     << ","
+     << dt
+     << ")\n";
+
+  CMPILCmdBuffer cmd_buffer(m_global_comm,m_global_rank);
+  CMPIBarrier barrier(m_global_comm);
+
+  //send cmd to slave
+  cmd_buffer.broadcast(CMD_ADD_SFIF);
+
+  AFieldMaster* new_fm =
+    new ScalarFluidInteractionFieldMaster(
+      &m_tml_global_comm,
+      fieldname,
+      filename,
+      savetype,
+      t_0,
+      t_end,
+      dt
+    );
+  m_save_fields.push_back(new_fm);
+
+  barrier.wait("CLatticeMaster::addScalarFluiedInteractionSaveField");
+
+  console.Debug() << "end CLatticeMaster::addScalarFluidInteractionSaveField()\n";
+}
+
+/*!
+  add a vector fluid interaction field to the list of fields to be saved
+
+  \param filename the name of the file the field is saved into
+  \param fieldname the name of the field
+  \param savetype output file format
+  \param t_0 first timestep to be saved
+  \param t_end last timestep to be saved
+  \param dt timesteps between saves
+*/
+void CLatticeMaster::addVectorFluidInteractionSaveField(
+  const string& filename,
+  const string& fieldname,
+  const string& savetype,
+  int t_0,
+  int t_end,
+  int dt
+)
+{
+  console.Debug()
+     << "CLatticeMaster::addVectorFluidInteractionSaveField("
+     << filename
+     << ","
+     << fieldname
+     << ","
+     << t_0
+     << ","
+     << t_end
+     << ","
+     << dt
+     << ")\n";
+
+  CMPILCmdBuffer cmd_buffer(m_global_comm,m_global_rank);
+  CMPIBarrier barrier(m_global_comm);
+
+  //send cmd to slave
+  cmd_buffer.broadcast(CMD_ADD_VFIF);
+
+  AFieldMaster* new_fm =
+    new VectorFluidInteractionFieldMaster(
+      &m_tml_global_comm,
+      fieldname,
+      filename,
+      savetype,
+      t_0,
+      t_end,
+      dt
+    );
+  m_save_fields.push_back(new_fm);
+
+  barrier.wait("CLatticeMaster::addVectorFluiedInteractionSaveField");
+
+  console.Debug() << "end CLatticeMaster::addVectorFluidInteractionSaveField()\n";
+}
+
+/****fluid contents: end****/
+
+
 void CLatticeMaster::init()
 {
   MPI_Comm_size(MPI_COMM_WORLD, &m_global_size);
@@ -151,7 +440,7 @@
 
   // -- local communicator = global comm - proc(0)
   // get global MPI_Group
-  MPI_Group mpi_global_group; 
+  MPI_Group mpi_global_group;
   MPI_Comm_group(MPI_COMM_WORLD,&mpi_global_group);
   // subtract id 0 from global group
   int id0=0;
@@ -163,12 +452,12 @@
   MPI_Group_free(&mpi_global_group);
 
   m_tml_global_comm.barrier();
-  	
+
   if(m_global_size!=ns+1){
     cerr << "wrong number of processes !! aborting" << endl << flush;
     exit(255);
     // should send abort command to workers
-  } 
+  }
   initializeConsole("console.out",1000); // 1st barrier -> sync time
 }
 
@@ -249,7 +538,7 @@
 
   //send command to slave
   cmd_buffer.broadcast(CMD_MAKELATTICE);
-  
+
   // send lattice parameters
 
   CLatticeParam clp(ptype, nrange, alpha, getProcessDims());
@@ -305,7 +594,7 @@
 
 /*!
  * Provides the initial minimum and maximum extents of all the particles read in from a geometry file.
- * 
+ *
  * \param initMinPt Minimum extent of particles inside domain.
  * \param initMaxPt Maximum extent of particles inside domain.
  */
@@ -328,10 +617,10 @@
 
   // check if domain has already been set
   if(m_geometry_is_initialized){
-	console.Warning() << "Spatial Domain has already been set and can not be changed - this call does nothing.\n"; 
-  } else { 	
-	// set geometry info according to parameters 
-	// no need to set circular bc - defaults to (0,0,0) anyway  
+	console.Warning() << "Spatial Domain has already been set and can not be changed - this call does nothing.\n";
+  } else {
+	// set geometry info according to parameters
+	// no need to set circular bc - defaults to (0,0,0) anyway
 	m_geo_info.setBBox(minBBoxPt,maxBBoxPt);
 	m_bbx_has_been_set=true;
 	// none of the boundaries circular -> simple initLattice
@@ -367,23 +656,23 @@
 
   // check if domain has already been set
   if(m_geometry_is_initialized){
-	console.Warning() << "Spatial Domain has already been set and can not be changed - this call does nothing.\n"; 
+	console.Warning() << "Spatial Domain has already been set and can not be changed - this call does nothing.\n";
   } else {
 	// check if we actually need the circular part or if the periodic dimensions are (false,false,false)
 	if((circDimVector[0]!=0) || (circDimVector[1]!=0) || (circDimVector[2]!=0)) { // we have at least one periodic dimension
-		// set geometry info according to parameters