diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 4d4e216c..e50a1aee 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -611,7 +611,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, //*************************************************************************************** double MorphDrain(DoubleArray &SignDist, signed char *id, - std::shared_ptr Dm, double VoidFraction) { + std::shared_ptr Dm, double VoidFraction, double InitialRadius) { // SignDist is the distance to the object that you want to constaing the morphological opening // VoidFraction is the the empty space where the object inst // id is a labeled map @@ -688,6 +688,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, double deltaR = 0.05; // amount to change the radius in voxel units double Rcrit_old = maxdistGlobal; double Rcrit_new = maxdistGlobal; + + if (InitialRadius < maxdistGlobal){ + Rcrit_old = InitialRadius; + Rcrit_new = InitialRadius; + } //if (argc>2){ // Rcrit_new = strtod(argv[2],NULL); // if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new); diff --git a/analysis/morphology.h b/analysis/morphology.h index 4ef2301c..681ee6dd 100644 --- a/analysis/morphology.h +++ b/analysis/morphology.h @@ -7,7 +7,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr Dm, double VoidFraction, signed char ErodeLabel, signed char ReplaceLabel); double MorphDrain(DoubleArray &SignDist, signed char *id, - std::shared_ptr Dm, double VoidFraction); + std::shared_ptr Dm, double VoidFraction, double InitialRadius); double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array &id, std::shared_ptr Dm, double TargetVol, double WallFactor); diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index a8e24273..8712c11e 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -53,6 +53,7 @@ int main(int argc, char **argv) auto WriteValues = domain_db->getVector( "WriteValues" ); SW = domain_db->getScalar("Sw"); auto READFILE = domain_db->getScalar( "Filename" ); + auto MORPH_RADIUS = domain_db->getWithDefault( "MorphRadius", 100000.0); // Generate the NWP configuration //if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit); @@ -122,7 +123,7 @@ int main(int argc, char **argv) comm.barrier(); // Run the morphological opening - MorphDrain(SignDist, id, Dm, SW); + MorphDrain(SignDist, id, Dm, SW, MORPH_RADIUS); // calculate distance to non-wetting fluid if (domain_db->keyExists( "HistoryLabels" )){