44 #include <boost/math/special_functions/fpclassify.hpp>
54 ProcessIsoContour::create,
55 "Extract an isocontour of fieldid variable and at "
56 "value fieldvalue, Optionally fieldstr can be "
57 "specified for a string defiition or smooth for "
65 "string of isocontour to be extracted");
67 "name for isocontour if fieldstr "
68 "specified, default is isocon");
71 "field id to extract");
74 "field value to extract");
77 "Globally condense contour to unique "
81 "Smooth isocontour (might require "
85 "Number of smoothing cycle, default = "
89 "Postive diffusion coefficient "
90 "(0 < lambda < 1), default = 0.5");
93 "Negative diffusion coefficient "
94 "(0 < mu < 1), default = 0.495");
97 "Remove contours with less than specified number of triangles."
98 "Only valid with GlobalCondense or Smooth options.");
110 int rank =
m_f->m_comm->GetRank();
116 cout <<
"Process Contour extraction..." << endl;
121 vector<IsoSharedPtr> iso;
123 if(
m_f->m_fieldPts.get())
129 if(
m_f->m_exp.size() == 0)
142 fieldid =
m_f->m_fieldPts->GetNFields();
148 string varstr =
"x y z";
149 vector<Array<OneD, const NekDouble> > interpfields;
151 for(
int i = 0; i <
m_f->m_fieldPts->GetDim(); ++i)
153 interpfields.push_back(
m_f->m_fieldPts->GetPts(i));
155 for(
int i = 0; i <
m_f->m_fieldPts->GetNFields(); ++i)
157 varstr +=
" " +
m_f->m_fieldPts->GetFieldName(i);
158 interpfields.push_back(
m_f->m_fieldPts->GetPts(i+3));
162 std::string str =
m_config[
"fieldstr"].as<
string>();
163 ExprId = strEval.DefineFunction(varstr.c_str(), str);
165 strEval.Evaluate(ExprId, interpfields, pts);
168 string fieldName =
m_config[
"fieldname"].as<
string>();
170 m_f->m_fieldPts->AddField(pts, fieldName);
174 ASSERTL0(
m_config[
"fieldid"].as<string>() !=
"NotSet",
"fieldid must be specified");
175 fieldid =
m_config[
"fieldid"].as<
int>();
178 ASSERTL0(
m_config[
"fieldvalue"].as<string>() !=
"NotSet",
"fieldvalue must be specified");
186 bool smoothing =
m_config[
"smooth"].m_beenSet;
187 bool globalcondense =
m_config[
"globalcondense"].m_beenSet;
190 int nfields =
m_f->m_fieldPts->GetNFields() +
m_f->m_fieldPts->GetDim();
193 g_iso->globalcondense(iso,
m_f->m_verbose);
197 iso.push_back(g_iso);
202 int niter =
m_config[
"smoothiter"].as<
int>();
205 for(
int i =0 ; i < iso.size(); ++i)
207 iso[i]->smooth(niter,lambda,-mu);
212 if((mincontour =
m_config[
"removesmallcontour"].as<int>()))
214 vector<IsoSharedPtr> new_iso;
218 cout <<
"Identifying separate regions [." << flush ;
220 for(
int i =0 ; i < iso.size(); ++i)
222 iso[i]->separate_regions(new_iso,mincontour,
m_f->m_verbose);
227 cout <<
"]" << endl << flush ;
245 ss << cpuTime <<
"s";
246 cout <<
"Process Isocontour CPU Time: " << setw(8) << left
265 if (((cx[0]-cx[1])==0.0)&&
266 ((cy[0]-cy[1])==0.0)&&
267 ((cz[0]-cz[1])==0.0))
269 if (((cx[2]-cx[3])==0.0)&&
270 ((cy[2]-cy[3])==0.0)&&
271 ((cz[2]-cz[3])==0.0))
350 ASSERTL0(
false,
"Error in 5-point triangulation in ThreeSimilar");
359 vector<IsoSharedPtr> returnval;
361 int coordim =
m_f->m_fieldPts->GetDim();
362 int nfields =
m_f->m_fieldPts->GetNFields() + coordim;
365 "This methods is currently only set up for 3D fields");
366 ASSERTL1(coordim + fieldid < nfields,
367 "field id is larger than number contained in FieldPts");
369 m_f->m_fieldPts->GetPts(fields);
373 int i, j, k, ii, jj, kk, r, s, n, counter, boolean;
376 for(i = 1; i < nfields; ++i)
378 intfields[i] = intfields[i-1] + 5;
384 vector<Array<OneD, int> > ptsConn;
385 m_f->m_fieldPts->GetConnectivity(ptsConn);
386 for(
int zone = 0; zone < ptsConn.size(); ++zone)
392 int nelmt = ptsConn[zone].num_elements()
397 for (n = 0, i = 0; i < nelmt; ++i)
400 if(!(((c[conn[i*4]] >val)&&(c[conn[i*4+1]]>val)&&
401 (c[conn[i*4+2]]>val)&&(c[conn[i*4+3]]>val))||
402 ((c[conn[i*4 ]]<val)&&(c[conn[i*4+1]]<val)&&
403 (c[conn[i*4+2]]<val)&&(c[conn[i*4+3]]<val))))
408 for (counter = 0, j=0; j<=2; j++)
410 for (k=j+1; k<=3; k++)
412 if (((c[conn[i*4+j]]>=val)&&
413 (val>=c[conn[i*4+k]]))||
414 ((c[conn[i*4+j]]<=val)&&
415 (val<=c[conn[i*4+k]])))
423 if(fabs(cj-ck) > 1e-12)
426 for(
int f = 0; f < nfields; ++f)
432 intfields[f][counter] =
433 fields[f][conn[4*i+j]] +
434 factor*(fields[f][conn[4*i+k]] -
435 fields[f][conn[4*i+j]]);
447 iso->resize_fields(3*n);
449 for(j = 0; j < 3; ++j)
451 iso->set_fields(3*(n-1)+j,intfields,j);
456 iso->resize_fields(3*n);
458 for(j = 0; j < 3; ++j)
460 iso->set_fields(3*(n-2)+j,intfields,j);
461 iso->set_fields(3*(n-1)+j,intfields,j+1);
466 iso->resize_fields(3*n);
469 for (ii=0;ii<=2;ii++)
471 for (jj=ii+1;jj<=3;jj++)
473 for (kk=jj+1;kk<=4;kk++)
475 if((((cx[ii]-cx[jj])==0.0)&&
476 ((cy[ii]-cy[jj])==0.0)&&
477 ((cz[ii]-cz[jj])==0.0))&&
478 (((cx[ii]-cx[kk])==0.0)&&
479 ((cy[ii]-cy[kk])==0.0)&&
480 ((cz[ii]-cz[kk])==0.0)))
485 iso->set_fields(3*(n-1) ,intfields,ii);
486 iso->set_fields(3*(n-1)+1,intfields,r);
487 iso->set_fields(3*(n-1)+2,intfields,s);
501 iso->set_fields(3*(n-1) ,intfields,0);
502 iso->set_fields(3*(n-1)+1,intfields,2);
503 iso->set_fields(3*(n-1)+2,intfields,r);
513 returnval.push_back(iso);
522 int nfields =
m_f->m_fieldPts->GetNFields() +
m_f->m_fieldPts->GetDim();
531 for(
int i =0; i < iso.size(); ++i)
533 npts += iso[i]->get_nvert();
542 for(
int i =0; i < iso.size(); ++i)
544 for(
int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
546 newfields[0][cnt] = iso[i]->get_x(j);
547 newfields[1][cnt] = iso[i]->get_y(j);
548 newfields[2][cnt] = iso[i]->get_z(j);
554 for(
int f = 0; f < nfields-3; ++f)
559 for(
int i =0; i < iso.size(); ++i)
561 for(
int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
563 newfields[f+3][cnt] = iso[i]->get_fields(f,j);
568 m_f->m_fieldPts->SetPts(newfields);
571 vector<Array<OneD, int> > ptsConn;
572 m_f->m_fieldPts->GetConnectivity(ptsConn);
575 for(
int i =0; i < iso.size(); ++i)
577 int ntris = iso[i]->get_ntris();
580 for(
int j = 0; j < 3*ntris; ++j)
582 conn[j] = cnt + iso[i]->get_vid(j);
584 ptsConn.push_back(conn);
585 cnt += iso[i]->get_nvert();
587 m_f->m_fieldPts->SetConnectivity(ptsConn);
594 "Assume input is from ePtsTriBlock");
597 int dim =
m_f->m_fieldPts->GetDim();
598 int nfields =
m_f->m_fieldPts->GetNFields() + dim;
600 m_f->m_fieldPts->GetPts(fieldpts);
601 vector<Array<OneD, int> > ptsConn;
602 m_f->m_fieldPts->GetConnectivity(ptsConn);
606 for(
int c = 0; c < ptsConn.size(); ++c)
612 nelmt = ptsConn[c].num_elements()/3;
614 iso->set_ntris(nelmt);
615 iso->resize_vid(3*nelmt);
619 for(
int i = 0; i < ptsConn[c].num_elements(); ++i)
621 int cid = ptsConn[c][i]-cnt;
623 nvert = max(cid,nvert);
627 iso->set_nvert(nvert);
628 iso->resize_fields(nvert);
631 for(
int i = 0; i < nvert; ++i)
633 iso->set_fields(i,fieldpts,i+cnt);
636 isovec.push_back(iso);
644 register int i,j,cnt;
646 vector<IsoVertex>
vert;
660 for(cnt =0, i = 0; i < 3; ++i)
665 for(
int f = 0; f <
m_fields.size(); ++f)
677 for(j = 0; j < 3; ++j)
683 pt =
find(vert.begin(),vert.end(),v);
686 m_vid[3*i+j] = pt[0].m_id;
692 for(
int f = 0; f <
m_fields.size(); ++f)
713 for(j = 3*i; j < 3*(m_ntris-1); ++j)
731 for(
int f = 0; f <
m_fields.size(); ++f)
738 m_x[i] = vert[i].m_x;
739 m_y[i] = vert[i].m_y;
740 m_z[i] = vert[i].m_z;
741 for(
int f = 0; f <
m_fields.size(); ++f)
743 m_fields[f][i] = vert[i].m_fields[f];
769 if((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2) < SQ_PNT_TOL)
791 for(i = 0; i < niso; ++i)
795 m_ntris += iso[i]->m_ntris;
802 for(i = 0; i < niso; ++i)
806 m_nvert += iso[i]->m_nvert;
810 vector< vector<int> > isocon;
817 vector<Array<OneD, NekDouble> > sph(niso);
819 for(i = 0; i < niso; ++i)
824 rng[0] = rng[3] = iso[i]->m_x[0];
825 rng[1] = rng[4] = iso[i]->m_x[1];
826 rng[2] = rng[5] = iso[i]->m_x[2];
828 for(id1 = 1; id1 < iso[i]->m_nvert;++id1)
830 rng[0] = min(rng[0],iso[i]->
m_x[i]);
831 rng[1] = min(rng[1],iso[i]->
m_y[i]);
832 rng[2] = min(rng[2],iso[i]->
m_z[i]);
834 rng[3] = max(rng[3],iso[i]->
m_x[i]);
835 rng[4] = max(rng[4],iso[i]->
m_y[i]);
836 rng[5] = max(rng[5],iso[i]->
m_z[i]);
840 sph[i][0] = (rng[3]+rng[0])/2.0;
841 sph[i][1] = (rng[4]+rng[1])/2.0;
842 sph[i][2] = (rng[5]+rng[2])/2.0;
845 sph[i][3] = sqrt((rng[3]-rng[0])*(rng[3]-rng[0]) +
846 (rng[4]-rng[1])*(rng[4]-rng[1]) +
847 (rng[5]-rng[2])*(rng[5]-rng[2]));
850 for(i = 0; i < niso; ++i)
852 for(j = i; j < niso; ++j)
854 NekDouble diff=sqrt((sph[i][0]-sph[j][0])*(sph[i][0]-sph[j][0])+
855 (sph[i][1]-sph[j][1])*(sph[i][1]-sph[j][1])+
856 (sph[i][2]-sph[j][2])*(sph[i][2]-sph[j][2]));
859 if(diff < sph[i][3] + sph[j][3])
861 isocon[i].push_back(j);
869 for(i = 0; i < niso; ++i)
877 for(i = 0; i < niso; ++i)
879 totiso += isocon[i].size();
885 cout <<
"Progress Bar totiso: " << totiso << endl;
887 for(i = 0; i < niso; ++i)
889 for(n = 0; n < isocon[i].size(); ++n, ++cnt)
892 if(verbose && totiso >= 40)
897 int con = isocon[i][n];
898 for(id1 = 0; id1 < iso[i]->m_nvert; ++id1)
901 if(verbose && totiso < 40)
911 for(id2 = start; id2 < iso[con]->m_nvert; ++id2)
914 if((vidmap[con][id2] == -1)||(vidmap[i][id1] == -1))
917 iso[i]->
m_z[id1], iso[con]->
m_x[id2],
918 iso[con]->
m_y[id2],iso[con]->
m_z[id2]))
920 if((vidmap[i][id1] == -1) &&
921 (vidmap[con][id2] != -1))
923 vidmap[i][id1] = vidmap[con][id2];
925 else if((vidmap[con][id2] == -1) &&
926 (vidmap[i][id1] != -1))
928 vidmap[con][id2] = vidmap[i][id1];
930 else if((vidmap[con][id2] == -1) &&
931 (vidmap[i][id1] == -1))
933 vidmap[i][id1] = vidmap[con][id2] = nvert++;
941 for(id1 = 0; id1 < iso[i]->m_nvert;++id1)
943 if(vidmap[i][id1] == -1)
945 vidmap[i][id1] = nvert++;
953 for(n = 0; n < niso; ++n)
955 for(i = 0; i < iso[n]->m_ntris; ++i,nelmt++)
959 m_vid[3*nelmt+j] = vidmap[n][iso[n]->m_vid[3*i+j]];
971 for(i = 0; i < iso[0]->m_fields.size(); ++i)
977 for(n = 0; n < niso; ++n)
979 for(i = 0; i < iso[n]->m_nvert; ++i)
981 m_x[vidmap[n][i]] = iso[n]->m_x[i];
982 m_y[vidmap[n][i]] = iso[n]->m_y[i];
983 m_z[vidmap[n][i]] = iso[n]->m_z[i];
985 for(j = 0; j <
m_fields.size(); ++j)
987 m_fields[j][vidmap[n][i]] = iso[n]->m_fields[j][i];
1000 vector< vector<int> > adj,vertcon;
1008 for(j = 0; j < 3; ++j)
1010 vertcon[
m_vid[3*i+j]].push_back(i);
1019 for(ipt = vertcon[i].begin(); ipt != vertcon[i].end(); ++ipt)
1021 for(j = 0; j < 3; ++j)
1024 if(
m_vid[3*(*ipt)+j] != i)
1027 for(iad = adj[i].begin(); iad != adj[i].end();++iad)
1029 if(*iad == (
m_vid[3*(*ipt)+j]))
break;
1032 if(iad == adj[i].end())
1034 adj[i].push_back(
m_vid[3*(*ipt)+j]);
1046 for (iter=0;iter<n_iter;iter++)
1053 del_v[0] = del_v[1] = del_v[2] = 0.0;
1055 for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
1057 del_v[0] = del_v[0] + (
m_x[*iad]-
m_x[i])*w;
1058 del_v[1] = del_v[1] + (
m_y[*iad]-
m_y[i])*w;
1059 del_v[2] = del_v[2] + (
m_z[*iad]-
m_z[i])*w;
1062 xtemp[i] =
m_x[i] + del_v[0] * lambda;
1063 ytemp[i] =
m_y[i] + del_v[1] * lambda;
1064 ztemp[i] =
m_z[i] + del_v[2] * lambda;
1072 del_v[0] = del_v[1] = del_v[2] = 0;
1074 for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
1076 del_v[0] = del_v[0] + (
m_x[*iad]-
m_x[i])*w;
1077 del_v[1] = del_v[1] + (
m_y[*iad]-
m_y[i])*w;
1078 del_v[2] = del_v[2] + (
m_z[*iad]-
m_z[i])*w;
1081 m_x[i] = xtemp[i] + del_v[0] * mu;
1082 m_y[i] = ytemp[i] + del_v[1] * mu;
1083 m_z[i] = ztemp[i] + del_v[2] * mu;
1102 for(j = 0; j < 3; ++j)
1104 vertcon[
m_vid[3*i+j]].push_back(i);
1121 if(vidregion[k] == -1)
1123 vidregion[k] = ++nregions;
1126 for(ipt = vertcon[k].begin(); ipt != vertcon[k].end(); ++ipt)
1128 for(i = 0; i < 3; ++i)
1130 if(vidregion[
id =
m_vid[3*(ipt[0])+i]] == -1)
1132 tocheck.push_back(
id);
1133 vidregion[id] = nregions;
1140 while(tocheck.size())
1142 cid = tocheck.begin();
1143 while(cid != tocheck.end())
1147 for(ipt = vertcon[*cid].begin(); ipt != vertcon[*cid].end(); ++ipt)
1149 for(i = 0; i < 3; ++i)
1151 if(vidregion[
id =
m_vid[3*(ipt[0])+i]] == -1)
1153 tocheck.push_back(
id);
1154 vidregion[id] = nregions;
1160 tocheck.pop_front();
1175 nvert[vidregion[i]] +=1;
1181 nelmt[vidregion[
m_vid[3*i]]] +=1;
1186 for(
int n = 0; n < nregions; ++n)
1188 if(nelmt[n] > minsize)
1193 iso->m_ntris = nelmt[n];
1196 iso->m_nvert = nvert[n];
1197 iso->m_x.resize(nvert[n]);
1198 iso->m_y.resize(nvert[n]);
1199 iso->m_z.resize(nvert[n]);
1201 iso->m_fields.resize(nfields);
1202 for(i = 0; i < nfields; ++i)
1204 iso->m_fields[i].resize(nvert[n]);
1213 if(vidregion[i] == n)
1222 if(vidregion[
m_vid[3*i]] == n)
1224 for(j = 0; j < 3; ++j)
1226 iso->m_vid[3*cnt+j] = vidmap[
m_vid[3*i+j]];
1228 iso->m_x[vidmap[m_vid[3*i+j]]] =
m_x[m_vid[3*i+j]];
1229 iso->m_y[vidmap[m_vid[3*i+j]]] =
m_y[m_vid[3*i+j]];
1230 iso->m_z[vidmap[m_vid[3*i+j]]] =
m_z[m_vid[3*i+j]];
1232 for(k = 0; k < nfields; ++k)
1234 iso->m_fields[k][vidmap[m_vid[3*i+j]]] =
m_fields[k][m_vid[3*i+j]];
1241 WARNINGL0(cnt == nelmt[n],
"Number of elements do not match");
1242 sep_iso.push_back(iso);
#define ASSERTL0(condition, msg)
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void SetupEquiSpacedField(void)
map< string, ConfigOption > m_config
List of configuration values.
virtual ~ProcessIsoContour()
FieldSharedPtr m_f
Field object.
void separate_regions(vector< boost::shared_ptr< Iso > > &iso, int minsize, bool verbose)
vector< NekDouble > m_fields
void smooth(int n_iter, NekDouble lambda, NekDouble mu)
boost::shared_ptr< Iso > IsoSharedPtr
#define WARNINGL0(condition, msg)
bool operator!=(const IsoVertex &x, const IsoVertex &y)
void SetupIsoFromFieldPts(vector< IsoSharedPtr > &isovec)
void TwoPairs(Array< OneD, NekDouble > &cx, Array< OneD, NekDouble > &cy, Array< OneD, NekDouble > &cz, int &pr)
boost::shared_ptr< Field > FieldSharedPtr
vector< vector< NekDouble > > m_fields
This class defines evaluator of analytic (symbolic) mathematical expressions. Expressions are allowed...
Represents a command-line configuration option.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool operator==(TriFaceIDs const &p1, TriFaceIDs const &p2)
vector< IsoSharedPtr > ExtractContour(const int fieldid, const NekDouble val)
void globalcondense(vector< boost::shared_ptr< Iso > > &iso, bool verbose)
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
void ResetFieldPts(vector< IsoSharedPtr > &iso)
bool same(NekDouble x1, NekDouble y1, NekDouble z1, NekDouble x2, NekDouble y2, NekDouble z2)
void ThreeSimilar(const int i, const int j, const int k, int &pr, int &ps)
void Zero(int n, T *x, const int incx)
Zero vector.
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
This processing module interpolates one field to another.
void PrintProgressbar(const int position, const int goal, const string message)
Prints a progressbar.
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.