2016-11-07 17:48:25 -06:00
|
|
|
/*
|
|
|
|
* Pre-processor to generate signed distance function from segmented data
|
|
|
|
* segmented data should be stored in a raw binary file as 1-byte integer (type char)
|
|
|
|
* will output distance functions for phases
|
|
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
|
|
#include <stdlib.h>
|
|
|
|
#include <math.h>
|
|
|
|
#include <iostream>
|
|
|
|
#include <fstream>
|
|
|
|
#include <sstream>
|
|
|
|
#include "common/Array.h"
|
|
|
|
#include "common/Domain.h"
|
2018-11-02 11:53:09 -05:00
|
|
|
#include "analysis/distance.h"
|
2018-12-03 15:06:38 -06:00
|
|
|
#include "analysis/morphology.h"
|
2016-11-07 17:48:25 -06:00
|
|
|
|
|
|
|
//*************************************************************************
|
|
|
|
// Morpohologica pre-processor
|
|
|
|
// Initialize phase distribution using morphological approach
|
|
|
|
// Signed distance function is used to determine fluid configuration
|
|
|
|
//*************************************************************************
|
|
|
|
|
|
|
|
int main(int argc, char **argv)
|
|
|
|
{
|
|
|
|
// Initialize MPI
|
2021-01-05 17:43:44 -06:00
|
|
|
Utilities::startup( argc, argv );
|
2020-01-28 07:51:32 -06:00
|
|
|
Utilities::MPI comm( MPI_COMM_WORLD );
|
2021-01-05 17:43:44 -06:00
|
|
|
int rank = comm.getRank();
|
2018-06-18 12:11:18 -05:00
|
|
|
{
|
|
|
|
//.......................................................................
|
|
|
|
// Reading the domain information file
|
|
|
|
//.......................................................................
|
|
|
|
char LocalRankFilename[40];
|
|
|
|
|
|
|
|
string filename;
|
2019-10-03 10:55:07 -05:00
|
|
|
double SW,Rcrit_new;
|
2018-06-19 06:12:03 -05:00
|
|
|
if (argc > 1){
|
2018-06-18 12:11:18 -05:00
|
|
|
filename=argv[1];
|
|
|
|
Rcrit_new=0.f;
|
2018-06-19 06:12:03 -05:00
|
|
|
//SW=strtod(argv[2],NULL);
|
2016-11-07 17:48:25 -06:00
|
|
|
}
|
2018-06-18 12:11:18 -05:00
|
|
|
else ERROR("No input database provided\n");
|
2020-01-22 11:01:29 -06:00
|
|
|
NULL_USE( Rcrit_new );
|
2018-06-18 12:11:18 -05:00
|
|
|
// read the input database
|
|
|
|
auto db = std::make_shared<Database>( filename );
|
|
|
|
auto domain_db = db->getDatabase( "Domain" );
|
|
|
|
|
|
|
|
// Read domain parameters
|
2019-10-03 11:14:32 -05:00
|
|
|
auto READFILE = domain_db->getScalar<std::string>( "Filename" );
|
2018-06-18 12:11:18 -05:00
|
|
|
auto size = domain_db->getVector<int>( "n" );
|
|
|
|
auto nproc = domain_db->getVector<int>( "nproc" );
|
2019-03-29 06:26:07 -05:00
|
|
|
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
|
|
|
|
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
|
2018-06-19 06:12:03 -05:00
|
|
|
SW = domain_db->getScalar<double>("Sw");
|
2019-04-26 16:36:22 -05:00
|
|
|
signed char ErodeLabel=2;
|
|
|
|
signed char OpenLabel=1;
|
|
|
|
if (domain_db->keyExists( "OpenLabel" )){
|
2019-04-26 16:46:46 -05:00
|
|
|
OpenLabel = domain_db->getScalar<int>("OpenLabel");
|
2019-04-26 16:36:22 -05:00
|
|
|
}
|
|
|
|
if (domain_db->keyExists( "ErodeLabel" )){
|
2019-04-26 16:46:46 -05:00
|
|
|
ErodeLabel = domain_db->getScalar<int>("ErodeLabel");
|
2019-04-26 16:36:22 -05:00
|
|
|
}
|
2018-12-03 15:06:38 -06:00
|
|
|
// Generate the NWP configuration
|
|
|
|
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
|
|
|
|
if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
|
|
|
|
// GenerateResidual(id,nx,ny,nz,Saturation);
|
|
|
|
|
2020-01-22 11:01:29 -06:00
|
|
|
int nx = size[0];
|
|
|
|
int ny = size[1];
|
|
|
|
int nz = size[2];
|
2018-06-18 12:11:18 -05:00
|
|
|
|
2020-01-22 11:01:29 -06:00
|
|
|
size_t N = (nx+2)*(ny+2)*(nz+2);
|
2018-06-18 12:11:18 -05:00
|
|
|
|
|
|
|
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
|
2019-10-03 10:55:07 -05:00
|
|
|
std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
|
2018-06-19 12:23:10 -05:00
|
|
|
// std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
|
2020-01-22 11:01:29 -06:00
|
|
|
for (size_t n=0; n<N; n++) Dm->id[n]=1;
|
2018-06-18 12:11:18 -05:00
|
|
|
Dm->CommInit();
|
2018-06-19 12:23:10 -05:00
|
|
|
|
2019-03-29 06:26:07 -05:00
|
|
|
signed char *id;
|
|
|
|
id = new signed char [N];
|
2019-10-03 10:55:07 -05:00
|
|
|
Mask->Decomp(READFILE);
|
|
|
|
Mask->CommInit();
|
2022-01-12 13:31:12 -06:00
|
|
|
|
2019-10-03 10:55:07 -05:00
|
|
|
// Generate the NWP configuration
|
|
|
|
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
|
|
|
|
if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
|
|
|
|
// GenerateResidual(id,nx,ny,nz,Saturation);
|
2018-11-02 11:50:10 -05:00
|
|
|
|
2018-06-19 08:55:09 -05:00
|
|
|
nx+=2; ny+=2; nz+=2;
|
2018-11-02 11:50:10 -05:00
|
|
|
// Generate the signed distance map
|
|
|
|
// Initialize the domain and communication
|
|
|
|
Array<char> id_solid(nx,ny,nz);
|
2018-06-18 12:11:18 -05:00
|
|
|
DoubleArray SignDist(nx,ny,nz);
|
2018-11-02 11:50:10 -05:00
|
|
|
|
|
|
|
// Solve for the position of the solid phase
|
|
|
|
for (int k=0;k<nz;k++){
|
|
|
|
for (int j=0;j<ny;j++){
|
|
|
|
for (int i=0;i<nx;i++){
|
|
|
|
int n = k*nx*ny+j*nx+i;
|
2019-10-03 11:19:07 -05:00
|
|
|
id[n] = Mask->id[n];
|
2018-11-02 11:50:10 -05:00
|
|
|
// Initialize the solid phase
|
2019-10-03 11:19:07 -05:00
|
|
|
if (Mask->id[n] > 0){
|
|
|
|
id_solid(i,j,k) = 1;
|
|
|
|
}
|
|
|
|
else
|
|
|
|
id_solid(i,j,k) = 0;
|
2018-11-02 11:50:10 -05:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
// Initialize the signed distance function
|
|
|
|
for (int k=0;k<nz;k++){
|
|
|
|
for (int j=0;j<ny;j++){
|
|
|
|
for (int i=0;i<nx;i++){
|
|
|
|
// Initialize distance to +/- 1
|
|
|
|
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
|
2018-11-02 11:57:03 -05:00
|
|
|
CalcDist(SignDist,id_solid,*Dm);
|
2018-06-18 12:11:18 -05:00
|
|
|
|
2021-01-05 17:43:44 -06:00
|
|
|
comm.barrier();
|
2018-06-18 12:11:18 -05:00
|
|
|
|
2018-12-03 15:06:38 -06:00
|
|
|
// Run the morphological opening
|
2019-04-26 16:36:22 -05:00
|
|
|
MorphOpen(SignDist, id, Dm, SW, ErodeLabel, OpenLabel);
|
2018-12-09 07:05:23 -06:00
|
|
|
|
|
|
|
// calculate distance to non-wetting fluid
|
|
|
|
if (domain_db->keyExists( "HistoryLabels" )){
|
|
|
|
if (rank==0) printf("Relabel solid components that touch fluid 1 \n");
|
2019-03-29 06:26:07 -05:00
|
|
|
auto LabelList = domain_db->getVector<int>( "ComponentLabels" );
|
|
|
|
auto HistoryLabels = domain_db->getVector<int>( "HistoryLabels" );
|
2018-12-09 07:05:23 -06:00
|
|
|
size_t NLABELS=LabelList.size();
|
2018-12-10 07:45:47 -06:00
|
|
|
if (rank==0){
|
|
|
|
for (unsigned int idx=0; idx < NLABELS; idx++){
|
2019-03-29 06:26:07 -05:00
|
|
|
signed char VALUE = LabelList[idx];
|
|
|
|
signed char NEWVAL = HistoryLabels[idx];
|
2018-12-10 07:45:47 -06:00
|
|
|
printf(" Relabel component %d as %d \n", VALUE, NEWVAL);
|
|
|
|
}
|
|
|
|
}
|
2018-12-09 07:05:23 -06:00
|
|
|
for (int k=0;k<nz;k++){
|
|
|
|
for (int j=0;j<ny;j++){
|
|
|
|
for (int i=0;i<nx;i++){
|
|
|
|
int n = k*nx*ny+j*nx+i;
|
|
|
|
// Initialize the solid phase
|
2018-12-09 17:01:59 -06:00
|
|
|
if (id[n] == 1) id_solid(i,j,k) = 0;
|
|
|
|
else id_solid(i,j,k) = 1;
|
2018-12-09 07:05:23 -06:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
// Initialize the signed distance function
|
|
|
|
for (int k=0;k<nz;k++){
|
|
|
|
for (int j=0;j<ny;j++){
|
|
|
|
for (int i=0;i<nx;i++){
|
|
|
|
// Initialize distance to +/- 1
|
|
|
|
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
CalcDist(SignDist,id_solid,*Dm);
|
|
|
|
// re-label IDs near the non-wetting fluid
|
|
|
|
for (int k=0;k<nz;k++){
|
|
|
|
for (int j=0;j<ny;j++){
|
|
|
|
for (int i=0;i<nx;i++){
|
|
|
|
int n = k*nx*ny+j*nx+i;
|
2018-12-10 08:16:37 -06:00
|
|
|
signed char LOCVAL = id[n];
|
2018-12-09 07:05:23 -06:00
|
|
|
for (unsigned int idx=0; idx < NLABELS; idx++){
|
2019-03-29 06:26:07 -05:00
|
|
|
signed char VALUE=LabelList[idx];
|
|
|
|
signed char NEWVALUE=HistoryLabels[idx];
|
2018-12-10 08:02:01 -06:00
|
|
|
if (LOCVAL == VALUE){
|
2018-12-09 07:05:23 -06:00
|
|
|
idx = NLABELS;
|
2019-04-26 13:46:25 -05:00
|
|
|
if (SignDist(i,j,k) < 2.0){
|
2018-12-10 08:16:37 -06:00
|
|
|
id[n] = NEWVALUE;
|
|
|
|
}
|
2018-12-09 07:05:23 -06:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2018-06-19 12:23:10 -05:00
|
|
|
|
2018-06-19 08:26:22 -05:00
|
|
|
if (rank==0) printf("Writing ID file \n");
|
2018-06-19 09:35:08 -05:00
|
|
|
sprintf(LocalRankFilename,"ID.%05i",rank);
|
2018-11-02 11:57:03 -05:00
|
|
|
|
2018-11-02 12:00:26 -05:00
|
|
|
FILE *ID = fopen(LocalRankFilename,"wb");
|
2018-06-18 12:11:18 -05:00
|
|
|
fwrite(id,1,N,ID);
|
|
|
|
fclose(ID);
|
2019-10-03 10:55:07 -05:00
|
|
|
|
|
|
|
// write the geometry to a single file
|
|
|
|
for (int k=0;k<nz;k++){
|
|
|
|
for (int j=0;j<ny;j++){
|
|
|
|
for (int i=0;i<nx;i++){
|
|
|
|
int n = k*nx*ny+j*nx+i;
|
|
|
|
Mask->id[n] = id[n];
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2021-01-05 17:43:44 -06:00
|
|
|
comm.barrier();
|
2019-10-03 10:55:07 -05:00
|
|
|
|
2020-01-22 11:01:29 -06:00
|
|
|
auto filename2 = READFILE + ".morphopen.raw";
|
|
|
|
if (rank==0) printf("Writing file to: %s \n", filename2.data());
|
|
|
|
Mask->AggregateLabels(filename2);
|
2017-07-09 17:57:34 -05:00
|
|
|
}
|
2016-11-07 17:48:25 -06:00
|
|
|
|
2020-10-08 10:03:42 -05:00
|
|
|
Utilities::shutdown();
|
2016-11-07 17:48:25 -06:00
|
|
|
}
|