38 #include <boost/geometry.hpp>
39 #include <boost/geometry/geometries/point.hpp>
40 #include <boost/geometry/geometries/box.hpp>
41 #include <boost/geometry/index/rtree.hpp>
48 #include <boost/math/special_functions/fpclassify.hpp>
50 namespace bg = boost::geometry;
51 namespace bgi = boost::geometry::index;
63 ProcessIsoContour::create,
64 "Extract an isocontour of fieldid variable and at "
65 "value fieldvalue, Optionally fieldstr can be "
66 "specified for a string defiition or smooth for "
74 "string of isocontour to be extracted");
76 "name for isocontour if fieldstr "
77 "specified, default is isocon");
80 "field id to extract");
83 "field value to extract");
86 "Globally condense contour to unique "
90 "Smooth isocontour (might require "
94 "Number of smoothing cycle, default = "
98 "Postive diffusion coefficient "
99 "(0 < lambda < 1), default = 0.5");
102 "Negative diffusion coefficient "
103 "(0 < mu < 1), default = 0.505");
106 "Remove contours with less than specified number of triangles."
107 "Only valid with GlobalCondense or Smooth options.");
119 int rank =
m_f->m_comm->GetRank();
125 cout <<
"Process Contour extraction..." << endl;
130 vector<IsoSharedPtr> iso;
132 if(
m_f->m_fieldPts.get())
136 cout <<
"Process read iso from Field Pts" << endl;
143 if(
m_f->m_exp.size() == 0)
156 fieldid =
m_f->m_fieldPts->GetNFields();
162 string varstr =
"x y z";
163 vector<Array<OneD, const NekDouble> > interpfields;
165 for(
int i = 0; i <
m_f->m_fieldPts->GetDim(); ++i)
167 interpfields.push_back(
m_f->m_fieldPts->GetPts(i));
169 for(
int i = 0; i <
m_f->m_fieldPts->GetNFields(); ++i)
171 varstr +=
" " +
m_f->m_fieldPts->GetFieldName(i);
172 interpfields.push_back(
m_f->m_fieldPts->GetPts(i+3));
176 std::string str =
m_config[
"fieldstr"].as<
string>();
177 ExprId = strEval.DefineFunction(varstr.c_str(), str);
179 strEval.Evaluate(ExprId, interpfields, pts);
182 string fieldName =
m_config[
"fieldname"].as<
string>();
184 m_f->m_fieldPts->AddField(pts, fieldName);
188 ASSERTL0(
m_config[
"fieldid"].as<string>() !=
"NotSet",
"fieldid must be specified");
189 fieldid =
m_config[
"fieldid"].as<
int>();
192 ASSERTL0(
m_config[
"fieldvalue"].as<string>() !=
"NotSet",
"fieldvalue must be specified");
200 bool smoothing =
m_config[
"smooth"].m_beenSet;
201 bool globalcondense =
m_config[
"globalcondense"].m_beenSet;
206 cout <<
"Process global condense ..." << endl;
208 int nfields =
m_f->m_fieldPts->GetNFields() +
m_f->m_fieldPts->GetDim();
211 g_iso->GlobalCondense(iso,
m_f->m_verbose);
215 iso.push_back(g_iso);
226 cout <<
"Process Contour smoothing ..." << endl;
231 int niter =
m_config[
"smoothiter"].as<
int>();
234 for(
int i =0 ; i < iso.size(); ++i)
236 iso[i]->Smooth(niter,lambda,-mu);
247 ss << cpuTime <<
"s";
248 cout <<
"Process smooth CPU Time: " << setw(8) << left
256 if((mincontour =
m_config[
"removesmallcontour"].as<int>()))
258 vector<IsoSharedPtr> new_iso;
262 cout <<
"Identifying separate regions [." << flush ;
264 for(
int i =0 ; i < iso.size(); ++i)
266 iso[i]->SeparateRegions(new_iso,mincontour,
m_f->m_verbose);
271 cout <<
"]" << endl << flush ;
289 ss << cpuTime <<
"s";
290 cout <<
"Process Isocontour CPU Time: " << setw(8) << left
309 if (((cx[0]-cx[1])==0.0)&&
310 ((cy[0]-cy[1])==0.0)&&
311 ((cz[0]-cz[1])==0.0))
313 if (((cx[2]-cx[3])==0.0)&&
314 ((cy[2]-cy[3])==0.0)&&
315 ((cz[2]-cz[3])==0.0))
394 ASSERTL0(
false,
"Error in 5-point triangulation in ThreeSimilar");
403 vector<IsoSharedPtr> returnval;
405 int coordim =
m_f->m_fieldPts->GetDim();
406 int nfields =
m_f->m_fieldPts->GetNFields() + coordim;
409 "This methods is currently only set up for 3D fields");
410 ASSERTL1(coordim + fieldid < nfields,
411 "field id is larger than number contained in FieldPts");
413 m_f->m_fieldPts->GetPts(fields);
417 int i, j, k, ii, jj, kk, r, s, n, counter, boolean;
420 for(i = 1; i < nfields; ++i)
422 intfields[i] = intfields[i-1] + 5;
428 vector<Array<OneD, int> > ptsConn;
429 m_f->m_fieldPts->GetConnectivity(ptsConn);
430 for(
int zone = 0; zone < ptsConn.size(); ++zone)
436 int nelmt = ptsConn[zone].num_elements()
441 for (n = 0, i = 0; i < nelmt; ++i)
444 if(!(((c[conn[i*4]] >val)&&(c[conn[i*4+1]]>val)&&
445 (c[conn[i*4+2]]>val)&&(c[conn[i*4+3]]>val))||
446 ((c[conn[i*4 ]]<val)&&(c[conn[i*4+1]]<val)&&
447 (c[conn[i*4+2]]<val)&&(c[conn[i*4+3]]<val))))
452 for (counter = 0, j=0; j<=2; j++)
454 for (k=j+1; k<=3; k++)
456 if (((c[conn[i*4+j]]>=val)&&
457 (val>=c[conn[i*4+k]]))||
458 ((c[conn[i*4+j]]<=val)&&
459 (val<=c[conn[i*4+k]])))
467 if(fabs(cj-ck) > 1e-12)
470 for(
int f = 0; f < nfields; ++f)
476 intfields[f][counter] =
477 fields[f][conn[4*i+j]] +
478 factor*(fields[f][conn[4*i+k]] -
479 fields[f][conn[4*i+j]]);
491 iso->ResizeFields(3*n);
493 for(j = 0; j < 3; ++j)
495 iso->SetFields(3*(n-1)+j,intfields,j);
500 iso->ResizeFields(3*n);
502 for(j = 0; j < 3; ++j)
504 iso->SetFields(3*(n-2)+j,intfields,j);
505 iso->SetFields(3*(n-1)+j,intfields,j+1);
510 iso->ResizeFields(3*n);
513 for (ii=0;ii<=2;ii++)
515 for (jj=ii+1;jj<=3;jj++)
517 for (kk=jj+1;kk<=4;kk++)
519 if((((cx[ii]-cx[jj])==0.0)&&
520 ((cy[ii]-cy[jj])==0.0)&&
521 ((cz[ii]-cz[jj])==0.0))&&
522 (((cx[ii]-cx[kk])==0.0)&&
523 ((cy[ii]-cy[kk])==0.0)&&
524 ((cz[ii]-cz[kk])==0.0)))
529 iso->SetFields(3*(n-1) ,intfields,ii);
530 iso->SetFields(3*(n-1)+1,intfields,r);
531 iso->SetFields(3*(n-1)+2,intfields,s);
545 iso->SetFields(3*(n-1) ,intfields,0);
546 iso->SetFields(3*(n-1)+1,intfields,2);
547 iso->SetFields(3*(n-1)+2,intfields,r);
560 returnval.push_back(iso);
570 int nfields =
m_f->m_fieldPts->GetNFields() +
m_f->m_fieldPts->GetDim();
579 for(
int i =0; i < iso.size(); ++i)
581 npts += iso[i]->GetNVert();
590 for(
int i =0; i < iso.size(); ++i)
592 for(
int j = 0; j < iso[i]->GetNVert(); ++j,++cnt)
594 newfields[0][cnt] = iso[i]->GetX(j);
595 newfields[1][cnt] = iso[i]->GetY(j);
596 newfields[2][cnt] = iso[i]->GetZ(j);
602 for(
int f = 0; f < nfields-3; ++f)
607 for(
int i =0; i < iso.size(); ++i)
609 for(
int j = 0; j < iso[i]->GetNVert(); ++j,++cnt)
611 newfields[f+3][cnt] = iso[i]->GetFields(f,j);
616 m_f->m_fieldPts->SetPts(newfields);
619 vector<Array<OneD, int> > ptsConn;
620 m_f->m_fieldPts->GetConnectivity(ptsConn);
623 for(
int i =0; i < iso.size(); ++i)
625 int ntris = iso[i]->GetNTris();
628 for(
int j = 0; j < 3*ntris; ++j)
630 conn[j] = cnt + iso[i]->GetVId(j);
632 ptsConn.push_back(conn);
633 cnt += iso[i]->GetNVert();
635 m_f->m_fieldPts->SetConnectivity(ptsConn);
642 "Assume input is from ePtsTriBlock");
645 int dim =
m_f->m_fieldPts->GetDim();
646 int nfields =
m_f->m_fieldPts->GetNFields() + dim;
648 m_f->m_fieldPts->GetPts(fieldpts);
649 vector<Array<OneD, int> > ptsConn;
650 m_f->m_fieldPts->GetConnectivity(ptsConn);
654 for(
int c = 0; c < ptsConn.size(); ++c)
660 nelmt = ptsConn[c].num_elements()/3;
662 iso->SetNTris(nelmt);
663 iso->ResizeVId(3*nelmt);
667 for(
int i = 0; i < ptsConn[c].num_elements(); ++i)
669 int cid = ptsConn[c][i]-cnt;
671 nvert = max(cid,nvert);
675 iso->SetNVert(nvert);
676 iso->ResizeFields(nvert);
679 for(
int i = 0; i < nvert; ++i)
681 iso->SetFields(i,fieldpts,i+cnt);
684 isovec.push_back(iso);
692 register int i,j,cnt;
694 vector<IsoVertex> vert;
708 for(cnt =0, i = 0; i < 3; ++i)
713 for(
int f = 0; f <
m_fields.size(); ++f)
725 for(j = 0; j < 3; ++j)
731 pt =
find(vert.begin(),vert.end(),v);
734 m_vid[3*i+j] = pt[0].m_id;
740 for(
int f = 0; f <
m_fields.size(); ++f)
761 for(j = 3*i; j < 3*(m_ntris-1); ++j)
779 for(
int f = 0; f <
m_fields.size(); ++f)
786 m_x[i] = vert[i].m_x;
787 m_y[i] = vert[i].m_y;
788 m_z[i] = vert[i].m_z;
789 for(
int f = 0; f <
m_fields.size(); ++f)
791 m_fields[f][i] = vert[i].m_fields[f];
799 typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> BPoint;
800 typedef std::pair<BPoint, unsigned int> PointPair;
814 for(i = 0; i < niso; ++i)
818 m_ntris += iso[i]->m_ntris;
825 for(i = 0; i < niso; ++i)
829 m_nvert += iso[i]->m_nvert;
833 vector< vector<int> > isocon;
835 vector<int> global_to_unique_map;
836 global_to_unique_map.resize(m_nvert);
838 vector< pair<int,int> > global_to_iso_map;
839 global_to_iso_map.resize(m_nvert);
843 std::vector<PointPair> inPoints;
844 for(i = 0; i < niso; ++i)
846 for(id1 = 0; id1 < iso[i]->m_nvert; ++id1)
848 inPoints.push_back(PointPair(BPoint( iso[i]->
m_x[id1],
850 iso[i]->
m_z[id1]), id2));
851 global_to_unique_map[id2]=id2;
852 global_to_iso_map[id2] = make_pair(i,id1);
860 cout <<
"Process building tree ..." << endl;
864 bgi::rtree<PointPair, bgi::rstar<16> > rtree;
865 rtree.insert(inPoints.begin(), inPoints.end());
868 int unique_index = 0;
875 "Nearest verts",prog);
878 BPoint queryPoint = inPoints[i].first;
883 if(global_to_unique_map[i] < unique_index)
890 std::vector<PointPair> result;
891 rtree.query(bgi::nearest(queryPoint, 10), std::back_inserter(result));
897 for(id1 = 0; id1 < result.size(); ++id1)
899 NekDouble dist = bg::distance(queryPoint, result[id1].first);
902 id2 = result[id1].second;
905 if(global_to_unique_map[id2] <unique_index)
907 new_index = global_to_unique_map[id2];
913 new_index = unique_index;
918 global_to_unique_map[i] = new_index;
919 for(it = samept.begin(); it != samept.end(); ++it)
921 global_to_unique_map[*it] = new_index;
931 m_nvert = unique_index;
936 for(n = 0; n < niso; ++n)
938 for(i = 0; i < iso[n]->m_ntris; ++i,nelmt++)
942 m_vid[3*nelmt+j] = global_to_unique_map[iso[n]->m_vid[3*i+j]+cnt];
945 cnt += iso[n]->m_nvert;
955 for(i = 0; i < iso[0]->m_fields.size(); ++i)
960 for(n = 0; n < global_to_unique_map.size(); ++n)
963 m_x[global_to_unique_map[n]] = bg::get<0>(inPoints[n].first);
964 m_y[global_to_unique_map[n]] = bg::get<1>(inPoints[n].first);
965 m_z[global_to_unique_map[n]] = bg::get<2>(inPoints[n].first);
967 int isoid = global_to_iso_map[n].first;
968 int ptid = global_to_iso_map[n].second;
969 for(j = 0; j <
m_fields.size(); ++j)
971 m_fields[j][global_to_unique_map[n]] = iso[isoid]->
998 vector< vector<int> > adj,vertcon;
999 vector< vector<NekDouble > > wght;
1007 for(j = 0; j < 3; ++j)
1009 vertcon[
m_vid[3*i+j]].push_back(i);
1020 for(ipt = vertcon[i].begin(); ipt != vertcon[i].end(); ++ipt)
1022 for(j = 0; j < 3; ++j)
1025 if(
m_vid[3*(*ipt)+j] != i)
1028 for(k = 0; k < adj[i].size(); ++k)
1030 if(adj[i][k] ==
m_vid[3*(*ipt)+j])
1036 if(k == adj[i].size())
1038 adj[i].push_back(
m_vid[3*(*ipt)+j]);
1049 for(j = 0; j < adj[i].size(); ++j)
1051 wght[i].push_back(w);
1061 for (iter=0;iter<n_iter;iter++)
1066 del_v[0] = del_v[1] = del_v[2] = 0.0;
1067 for(j = 0; j < adj[i].size(); ++j)
1069 del_v[0] = del_v[0] + (
m_x[adj[i][j]]-
m_x[i])*wght[i][j];
1070 del_v[1] = del_v[1] + (
m_y[adj[i][j]]-
m_y[i])*wght[i][j];
1071 del_v[2] = del_v[2] + (
m_z[adj[i][j]]-
m_z[i])*wght[i][j];
1074 xtemp[i] =
m_x[i] + del_v[0] * lambda;
1075 ytemp[i] =
m_y[i] + del_v[1] * lambda;
1076 ztemp[i] =
m_z[i] + del_v[2] * lambda;
1082 del_v[0] = del_v[1] = del_v[2] = 0;
1083 for(j = 0; j < adj[i].size(); ++j)
1085 del_v[0] = del_v[0] + (xtemp[adj[i][j]]-xtemp[i])*wght[i][j];
1086 del_v[1] = del_v[1] + (ytemp[adj[i][j]]-ytemp[i])*wght[i][j];
1087 del_v[2] = del_v[2] + (ztemp[adj[i][j]]-ztemp[i])*wght[i][j];
1090 m_x[i] = xtemp[i] + del_v[0] * mu;
1091 m_y[i] = ytemp[i] + del_v[1] * mu;
1092 m_z[i] = ztemp[i] + del_v[2] * mu;
1111 for(j = 0; j < 3; ++j)
1113 vertcon[
m_vid[3*i+j]].push_back(i);
1131 if(vidregion[k] == -1)
1133 vidregion[k] = ++nregions;
1136 for(ipt = vertcon[k].begin(); ipt != vertcon[k].end(); ++ipt)
1138 for(i = 0; i < 3; ++i)
1140 if(vidregion[
id =
m_vid[3*(ipt[0])+i]] == -1)
1142 tocheck.push_back(
id);
1143 vidregion[id] = nregions;
1150 while(tocheck.size())
1152 cid = tocheck.begin();
1153 while(cid != tocheck.end())
1157 for(ipt = vertcon[*cid].begin(); ipt != vertcon[*cid].end(); ++ipt)
1159 for(i = 0; i < 3; ++i)
1161 if(vidregion[
id =
m_vid[3*(ipt[0])+i]] == -1)
1163 tocheck.push_back(
id);
1164 vidregion[id] = nregions;
1170 tocheck.pop_front();
1185 nvert[vidregion[i]] +=1;
1191 nelmt[vidregion[
m_vid[3*i]]] +=1;
1196 for(
int n = 0; n < nregions; ++n)
1198 if(nelmt[n] > minsize)
1203 iso->m_ntris = nelmt[n];
1206 iso->m_nvert = nvert[n];
1207 iso->m_x.resize(nvert[n]);
1208 iso->m_y.resize(nvert[n]);
1209 iso->m_z.resize(nvert[n]);
1211 iso->m_fields.resize(nfields);
1212 for(i = 0; i < nfields; ++i)
1214 iso->m_fields[i].resize(nvert[n]);
1223 if(vidregion[i] == n)
1232 if(vidregion[
m_vid[3*i]] == n)
1234 for(j = 0; j < 3; ++j)
1236 iso->m_vid[3*cnt+j] = vidmap[
m_vid[3*i+j]];
1238 iso->m_x[vidmap[m_vid[3*i+j]]] =
m_x[m_vid[3*i+j]];
1239 iso->m_y[vidmap[m_vid[3*i+j]]] =
m_y[m_vid[3*i+j]];
1240 iso->m_z[vidmap[m_vid[3*i+j]]] =
m_z[m_vid[3*i+j]];
1242 for(k = 0; k < nfields; ++k)
1244 iso->m_fields[k][vidmap[m_vid[3*i+j]]] =
m_fields[k][m_vid[3*i+j]];
1251 WARNINGL0(cnt == nelmt[n],
"Number of elements do not match");
1252 sep_iso.push_back(iso);
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
bool operator==(TriFaceIDs const &p1, TriFaceIDs const &p2)
vector< NekDouble > m_fields
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Represents a command-line configuration option.
pair< ModuleType, string > ModuleKey
void ResetFieldPts(vector< IsoSharedPtr > &iso)
void ThreeSimilar(const int i, const int j, const int k, int &pr, int &ps)
void SetupIsoFromFieldPts(vector< IsoSharedPtr > &isovec)
vector< vector< NekDouble > > m_fields
void TwoPairs(Array< OneD, NekDouble > &cx, Array< OneD, NekDouble > &cy, Array< OneD, NekDouble > &cz, int &pr)
bool operator!=(const IsoVertex &x, const IsoVertex &y)
static const NekDouble kNekZeroTol
virtual ~ProcessIsoContour()
boost::shared_ptr< Iso > IsoSharedPtr
This processing module interpolates one field to another.
boost::shared_ptr< Field > FieldSharedPtr
#define WARNINGL0(condition, msg)
void GlobalCondense(vector< boost::shared_ptr< Iso > > &iso, bool verbose)
void SetupEquiSpacedField(void)
void Smooth(int n_iter, NekDouble lambda, NekDouble mu)
This class defines evaluator of analytic (symbolic) mathematical expressions. Expressions are allowed...
virtual void Process(po::variables_map &vm)
Write mesh to output file.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void SeparateRegions(vector< boost::shared_ptr< Iso > > &iso, int minsize, bool verbose)
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
vector< IsoSharedPtr > ExtractContour(const int fieldid, const NekDouble val)
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...
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.