Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessIsoContour.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessIsoContour.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Generate isocontours from field data.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 #include <string>
36 #include <iostream>
37 
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>
42 
43 #include "ProcessIsoContour.h"
44 
48 #include <boost/math/special_functions/fpclassify.hpp>
49 
50 namespace bg = boost::geometry;
51 namespace bgi = boost::geometry::index;
52 
53 using namespace std;
54 
55 namespace Nektar
56 {
57 namespace FieldUtils
58 {
59 
60 ModuleKey ProcessIsoContour::className =
62  ModuleKey(eProcessModule, "isocontour"),
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 "
67  "smoothing");
68 
69 ProcessIsoContour::ProcessIsoContour(FieldSharedPtr f) :
71 {
72 
73  m_config["fieldstr"] = ConfigOption(false, "NotSet",
74  "string of isocontour to be extracted");
75  m_config["fieldname"] = ConfigOption(false, "isocon",
76  "name for isocontour if fieldstr "
77  "specified, default is isocon");
78 
79  m_config["fieldid"] = ConfigOption(false, "NotSet",
80  "field id to extract");
81 
82  m_config["fieldvalue"] = ConfigOption(false, "NotSet",
83  "field value to extract");
84 
85  m_config["globalcondense"] = ConfigOption(true, "NotSet",
86  "Globally condense contour to unique "
87  "values");
88 
89  m_config["smooth"] = ConfigOption(true, "NotSet",
90  "Smooth isocontour (might require "
91  "globalcondense)");
92 
93  m_config["smoothiter"] = ConfigOption(false, "100",
94  "Number of smoothing cycle, default = "
95  "100");
96 
97  m_config["smoothposdiffusion"] = ConfigOption(false,"0.5",
98  "Postive diffusion coefficient "
99  "(0 < lambda < 1), default = 0.5");
100 
101  m_config["smoothnegdiffusion"] = ConfigOption(false,"0.505",
102  "Negative diffusion coefficient "
103  "(0 < mu < 1), default = 0.505");
104 
105  m_config["removesmallcontour"] = ConfigOption(false,"0",
106  "Remove contours with less than specified number of triangles."
107  "Only valid with GlobalCondense or Smooth options.");
108 
109 
110 }
111 
113 {
114 }
115 
116 void ProcessIsoContour::Process(po::variables_map &vm)
117 {
118  Timer timer;
119  int rank = m_f->m_comm->GetRank();
120 
121  if(m_f->m_verbose)
122  {
123  if(rank == 0)
124  {
125  cout << "Process Contour extraction..." << endl;
126  timer.Start();
127  }
128  }
129 
130  vector<IsoSharedPtr> iso;
131 
132  if(m_f->m_fieldPts.get()) // assume we have read .dat file to directly input dat file.
133  {
134  if(rank == 0)
135  {
136  cout << "Process read iso from Field Pts" << endl;
137  }
138 
140  }
141  else // extract isocontour from field
142  {
143  if(m_f->m_exp.size() == 0)
144  {
145  return;
146  }
147 
148  // extract all fields to equi-spaced
150 
151  int fieldid;
152  NekDouble value;
153 
154  if(m_config["fieldstr"].m_beenSet) //generate field of interest
155  {
156  fieldid = m_f->m_fieldPts->GetNFields();
157 
158  Array<OneD, NekDouble> pts(m_f->m_fieldPts->GetNpoints());
159 
160  // evaluate new function
162  string varstr = "x y z";
163  vector<Array<OneD, const NekDouble> > interpfields;
164 
165  for(int i = 0; i < m_f->m_fieldPts->GetDim(); ++i)
166  {
167  interpfields.push_back(m_f->m_fieldPts->GetPts(i));
168  }
169  for(int i = 0; i < m_f->m_fieldPts->GetNFields(); ++i)
170  {
171  varstr += " " + m_f->m_fieldPts->GetFieldName(i);
172  interpfields.push_back(m_f->m_fieldPts->GetPts(i+3));
173  }
174 
175  int ExprId = -1;
176  std::string str = m_config["fieldstr"].as<string>();
177  ExprId = strEval.DefineFunction(varstr.c_str(), str);
178 
179  strEval.Evaluate(ExprId, interpfields, pts);
180 
181  // set up field name if provided otherwise called "isocon" from default
182  string fieldName = m_config["fieldname"].as<string>();
183 
184  m_f->m_fieldPts->AddField(pts, fieldName);
185  }
186  else
187  {
188  ASSERTL0(m_config["fieldid"].as<string>() != "NotSet", "fieldid must be specified");
189  fieldid = m_config["fieldid"].as<int>();
190  }
191 
192  ASSERTL0(m_config["fieldvalue"].as<string>() != "NotSet", "fieldvalue must be specified");
193  value = m_config["fieldvalue"].as<NekDouble>();
194 
195  iso = ExtractContour(fieldid,value);
196 
197  }
198 
199  // Process isocontour
200  bool smoothing = m_config["smooth"].m_beenSet;
201  bool globalcondense = m_config["globalcondense"].m_beenSet;
202  if(globalcondense)
203  {
204  if(rank == 0)
205  {
206  cout << "Process global condense ..." << endl;
207  }
208  int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
210 
211  g_iso->GlobalCondense(iso,m_f->m_verbose);
212 
213 
214  iso.clear();
215  iso.push_back(g_iso);
216  }
217 
218  if(smoothing)
219  {
220  Timer timersm;
221 
222  if(m_f->m_verbose)
223  {
224  if(rank == 0)
225  {
226  cout << "Process Contour smoothing ..." << endl;
227  timersm.Start();
228  }
229  }
230 
231  int niter = m_config["smoothiter"].as<int>();
232  NekDouble lambda = m_config["smoothposdiffusion"].as<NekDouble>();
233  NekDouble mu = m_config["smoothnegdiffusion"].as<NekDouble>();
234  for(int i =0 ; i < iso.size(); ++i)
235  {
236  iso[i]->Smooth(niter,lambda,-mu);
237  }
238 
239  if(m_f->m_verbose)
240  {
241  if(rank == 0)
242  {
243  timersm.Stop();
244  NekDouble cpuTime = timersm.TimePerTest(1);
245 
246  stringstream ss;
247  ss << cpuTime << "s";
248  cout << "Process smooth CPU Time: " << setw(8) << left
249  << ss.str() << endl;
250  cpuTime = 0.0;
251  }
252  }
253  }
254 
255  int mincontour = 0;
256  if((mincontour = m_config["removesmallcontour"].as<int>()))
257  {
258  vector<IsoSharedPtr> new_iso;
259 
260  if(rank == 0)
261  {
262  cout << "Identifying separate regions [." << flush ;
263  }
264  for(int i =0 ; i < iso.size(); ++i)
265  {
266  iso[i]->SeparateRegions(new_iso,mincontour,m_f->m_verbose);
267  }
268 
269  if(rank == 0)
270  {
271  cout << "]" << endl << flush ;
272  }
273 
274  // reset iso to new_iso;
275  iso = new_iso;
276  }
277 
278  ResetFieldPts(iso);
279 
280 
281  if(m_f->m_verbose)
282  {
283  if(rank == 0)
284  {
285  timer.Stop();
286  NekDouble cpuTime = timer.TimePerTest(1);
287 
288  stringstream ss;
289  ss << cpuTime << "s";
290  cout << "Process Isocontour CPU Time: " << setw(8) << left
291  << ss.str() << endl;
292  cpuTime = 0.0;
293  }
294  }
295 }
296 
297 
298 /*********************/
299 /* Function TwoPairs */
300 /*********************/
301 
305  int &pr)
306 /* The logic of the following SWITCH statement may only be
307  understood with a "peusdo truth-table" */
308 {
309  if (((cx[0]-cx[1])==0.0)&&
310  ((cy[0]-cy[1])==0.0)&&
311  ((cz[0]-cz[1])==0.0))
312  {
313  if (((cx[2]-cx[3])==0.0)&&
314  ((cy[2]-cy[3])==0.0)&&
315  ((cz[2]-cz[3])==0.0))
316  {
317  pr=4;
318  }
319  else
320  {
321  pr=3;
322  }
323  }
324  else
325  {
326  pr=1;
327  }
328 }
329 
330 
331 /*************************/
332 /* Function ThreeSimilar */
333 /*************************/
334 void ThreeSimilar (const int i, const int j, const int k,
335  int &pr, int &ps)
336 /* The logic of the following SWITCH statement may be
337  understood with a "peusdo truth-table" */
338 {
339  switch (i + j + k)
340  {
341  case (3):
342  pr = 3;
343  ps = 4;
344  break;
345  case (4):
346  pr = 2;
347  ps = 4;
348  break;
349  case (5):
350  if (j == 1)
351  {
352  pr = 2;
353  ps = 3;
354  }
355  else
356  {
357  pr = 1;
358  ps = 4;
359  }
360  break;
361  case (6):
362  if (i == 0)
363  {
364  pr = 1;
365  ps = 3;
366  }
367  else
368  {
369  pr = 0;
370  ps = 4;
371  }
372  break;
373  case (7):
374  if (i == 0)
375  {
376  pr = 1;
377  ps = 2;
378  }
379  else
380  {
381  pr = 0;
382  ps = 3;
383  }
384  break;
385  case (8):
386  pr = 0;
387  ps = 2;
388  break;
389  case (9):
390  pr = 0;
391  ps = 1;
392  break;
393  default:
394  ASSERTL0(false,"Error in 5-point triangulation in ThreeSimilar");
395  break;
396  }
397 }
398 
400  const int fieldid,
401  const NekDouble val)
402 {
403  vector<IsoSharedPtr> returnval;
404 
405  int coordim = m_f->m_fieldPts->GetDim();
406  int nfields = m_f->m_fieldPts->GetNFields() + coordim;
407 
408  ASSERTL0(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock,
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);
414 
415  Array<OneD, NekDouble> c = fields[coordim + fieldid];
416 
417  int i, j, k, ii, jj, kk, r, s, n, counter, boolean;
418  Array<OneD, Array<OneD, NekDouble> > intfields(nfields);
419  intfields[0] = Array<OneD, NekDouble>(5*nfields);
420  for(i = 1; i < nfields; ++i)
421  {
422  intfields[i] = intfields[i-1] + 5;
423  }
424  Array<OneD, NekDouble> cx = intfields[0];
425  Array<OneD, NekDouble> cy = intfields[1];
426  Array<OneD, NekDouble> cz = intfields[2];
427 
428  vector<Array<OneD, int> > ptsConn;
429  m_f->m_fieldPts->GetConnectivity(ptsConn);
430  for(int zone = 0; zone < ptsConn.size(); ++zone)
431  {
432  IsoSharedPtr iso;
433 
434  iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-3);
435 
436  int nelmt = ptsConn[zone].num_elements()
437  /(coordim+1);
438 
439  Array<OneD, int> conn = ptsConn[zone];
440 
441  for (n = 0, i = 0; i < nelmt; ++i)
442  {
443  // check to see if val is between vertex values
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))))
448  {
449 
450  // loop over all edges and interpolate if
451  // contour is between vertex values
452  for (counter = 0, j=0; j<=2; j++)
453  {
454  for (k=j+1; k<=3; k++)
455  {
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]])))
460  {
461  // linear interpolation of fields
462  // (and coords).
463  NekDouble cj = c[conn[i*4+j]];
464  NekDouble ck = c[conn[i*4+k]];
465  NekDouble factor = (val-cj)/(ck-cj);
466 
467  if(fabs(cj-ck) > 1e-12)
468  {
469  // interpolate coordinates and fields
470  for(int f = 0; f < nfields; ++f)
471  {
472  if(counter == 5)
473  {
474  ASSERTL0(false,"Counter is 5");
475  }
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]]);
480  }
481  ++counter;
482  }
483  }
484  }
485  }
486 
487  switch(counter)
488  {
489  case 3:
490  n+=1;
491  iso->ResizeFields(3*n);
492 
493  for(j = 0; j < 3; ++j)
494  {
495  iso->SetFields(3*(n-1)+j,intfields,j);
496  }
497  break;
498  case 4:
499  n+=2;
500  iso->ResizeFields(3*n);
501 
502  for(j = 0; j < 3; ++j)
503  {
504  iso->SetFields(3*(n-2)+j,intfields,j);
505  iso->SetFields(3*(n-1)+j,intfields,j+1);
506  }
507  break;
508  case 5:
509  n+=1;
510  iso->ResizeFields(3*n);
511 
512  boolean=0;
513  for (ii=0;ii<=2;ii++)
514  {
515  for (jj=ii+1;jj<=3;jj++)
516  {
517  for (kk=jj+1;kk<=4;kk++)
518  {
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)))
525  {
526  boolean+=1;
527  ThreeSimilar (ii,jj,kk,r,s);
528 
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);
532  }
533  else
534  {
535  boolean+=0;
536  }
537  }
538  }
539  }
540 
541  if (boolean==0)
542  {
543  TwoPairs (cx,cy,cz,r);
544 
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);
548  }
549  break;
550  }
551  }
552  }
553 
554  if(n)
555  {
556  iso->SetNTris(n);
557 
558  // condense the information in this elemental extraction.
559  iso->Condense();
560  returnval.push_back(iso);
561  }
562  }
563 
564  return returnval;
565 }
566 
567 // reset m_fieldPts with values from iso;
568 void ProcessIsoContour::ResetFieldPts(vector<IsoSharedPtr> &iso)
569 {
570  int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
571 
572  // set output to triangle block.
573  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
574 
575  Array<OneD, Array<OneD, NekDouble> > newfields(nfields);
576 
577  // count up number of points
578  int npts = 0;
579  for(int i =0; i < iso.size(); ++i)
580  {
581  npts += iso[i]->GetNVert();
582  }
583 
584  // set up coordinate in new field
585  newfields[0] = Array<OneD, NekDouble>(npts);
586  newfields[1] = Array<OneD, NekDouble>(npts);
587  newfields[2] = Array<OneD, NekDouble>(npts);
588 
589  int cnt = 0;
590  for(int i =0; i < iso.size(); ++i)
591  {
592  for(int j = 0; j < iso[i]->GetNVert(); ++j,++cnt)
593  {
594  newfields[0][cnt] = iso[i]->GetX(j);
595  newfields[1][cnt] = iso[i]->GetY(j);
596  newfields[2][cnt] = iso[i]->GetZ(j);
597  }
598  }
599 
600 
601  // set up fields
602  for(int f = 0; f < nfields-3; ++f)
603  {
604  newfields[f+3] = Array<OneD, NekDouble>(npts);
605 
606  cnt = 0;
607  for(int i =0; i < iso.size(); ++i)
608  {
609  for(int j = 0; j < iso[i]->GetNVert(); ++j,++cnt)
610  {
611  newfields[f+3][cnt] = iso[i]->GetFields(f,j);
612  }
613  }
614  }
615 
616  m_f->m_fieldPts->SetPts(newfields);
617 
618  // set up connectivity data.
619  vector<Array<OneD, int> > ptsConn;
620  m_f->m_fieldPts->GetConnectivity(ptsConn);
621  cnt = 0;
622  ptsConn.clear();
623  for(int i =0; i < iso.size(); ++i)
624  {
625  int ntris = iso[i]->GetNTris();
626  Array<OneD, int> conn(ntris*3);
627 
628  for(int j = 0; j < 3*ntris; ++j)
629  {
630  conn[j] = cnt + iso[i]->GetVId(j);
631  }
632  ptsConn.push_back(conn);
633  cnt += iso[i]->GetNVert();
634  }
635  m_f->m_fieldPts->SetConnectivity(ptsConn);
636 }
637 
638 // reset m_fieldPts with values from iso;
639 void ProcessIsoContour::SetupIsoFromFieldPts(vector<IsoSharedPtr> &isovec)
640 {
641  ASSERTL0(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock,
642  "Assume input is from ePtsTriBlock");
643 
644  // get information from PtsField
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);
651 
652 
653  int cnt = 0;
654  for(int c = 0; c < ptsConn.size(); ++c)
655  {
656  // set up single iso with all the information from PtsField
658 
659  int nelmt = 0;
660  nelmt = ptsConn[c].num_elements()/3;
661 
662  iso->SetNTris(nelmt);
663  iso->ResizeVId(3*nelmt);
664 
665  // fill in connectivity values.
666  int nvert = 0;
667  for(int i = 0; i < ptsConn[c].num_elements(); ++i)
668  {
669  int cid = ptsConn[c][i]-cnt;
670  iso->SetVId(i,cid);
671  nvert = max(cid,nvert);
672  }
673  nvert++;
674 
675  iso->SetNVert(nvert);
676  iso->ResizeFields(nvert);
677 
678  // fill in points values (including coordinates)
679  for(int i = 0; i < nvert; ++i)
680  {
681  iso->SetFields(i,fieldpts,i+cnt);
682  }
683  cnt += nvert;
684  isovec.push_back(iso);
685  }
686 
687 }
688 
689 
690 void Iso::Condense(void)
691 {
692  register int i,j,cnt;
693  IsoVertex v;
694  vector<IsoVertex> vert;
696 
697  if(!m_ntris) return;
698 
699  if(m_condensed) return;
700  m_condensed = true;
701 
702  vert.reserve(m_ntris/6);
703 
705 
706  // fill first 3 points and initialise fields
707  v.m_fields.resize(m_fields.size());
708  for(cnt =0, i = 0; i < 3; ++i)
709  {
710  v.m_x = m_x[i];
711  v.m_y = m_y[i];
712  v.m_z = m_z[i];
713  for(int f = 0; f < m_fields.size(); ++f)
714  {
715  v.m_fields[f] = m_fields[f][i];
716  }
717  v.m_id = cnt;
718  vert.push_back(v);
719  m_vid[i] = v.m_id;
720  ++cnt;
721  }
722 
723  for(i = 1; i < m_ntris; ++i)
724  {
725  for(j = 0; j < 3; ++j)
726  {
727  v.m_x = m_x[3*i+j];
728  v.m_y = m_y[3*i+j];
729  v.m_z = m_z[3*i+j];
730 
731  pt = find(vert.begin(),vert.end(),v);
732  if(pt != vert.end())
733  {
734  m_vid[3*i+j] = pt[0].m_id;
735  }
736  else
737  {
738  v.m_id = cnt;
739 
740  for(int f = 0; f < m_fields.size(); ++f)
741  {
742  v.m_fields[f] = m_fields[f][3*i+j];
743  }
744 
745  vert.push_back(v);
746 
747  m_vid[3*i+j] = v.m_id;
748  ++cnt;
749  }
750  }
751  }
752 
753  // remove elements with multiple vertices
754  for(i = 0,cnt=0; i < m_ntris;)
755  {
756  if((m_vid[3*i] ==m_vid[3*i+1])||
757  (m_vid[3*i] ==m_vid[3*i+2])||
758  (m_vid[3*i+1]==m_vid[3*i+2]))
759  {
760  cnt++;
761  for(j = 3*i; j < 3*(m_ntris-1); ++j)
762  {
763  m_vid[j] = m_vid[j+3];
764  }
765  m_ntris--;
766  }
767  else
768  {
769  ++i;
770  }
771  }
772 
773  m_nvert = vert.size();
774 
775  m_x.resize(m_nvert);
776  m_y.resize(m_nvert);
777  m_z.resize(m_nvert);
778 
779  for(int f = 0; f < m_fields.size(); ++f)
780  {
781  m_fields[f].resize(m_nvert);
782  }
783 
784  for(i = 0; i < m_nvert; ++i)
785  {
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)
790  {
791  m_fields[f][i] = vert[i].m_fields[f];
792  }
793  }
794 }
795 
796 
797 void Iso::GlobalCondense(vector<IsoSharedPtr> &iso, bool verbose)
798 {
799  typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> BPoint;
800  typedef std::pair<BPoint, unsigned int> PointPair;
801 
802  int i,j,n;
803  int nelmt;
804  int niso=iso.size();
805  int id1,id2;
807 
808  if(m_condensed) return;
809  m_condensed = true;
810 
811  vidmap = Array<OneD, Array<OneD, int> > (niso);
812 
813  m_ntris = 0;
814  for(i = 0; i < niso; ++i)
815  {
816  if(iso[i]->m_ntris)
817  {
818  m_ntris += iso[i]->m_ntris;
819  }
820  }
821 
823 
824  m_nvert = 0;
825  for(i = 0; i < niso; ++i)
826  {
827  if(iso[i]->m_ntris)
828  {
829  m_nvert += iso[i]->m_nvert;
830  }
831  }
832 
833  vector< vector<int> > isocon;
834  isocon.resize(niso);
835  vector<int> global_to_unique_map;
836  global_to_unique_map.resize(m_nvert);
837 
838  vector< pair<int,int> > global_to_iso_map;
839  global_to_iso_map.resize(m_nvert);
840 
841  //Fill vertex array into rtree format
842  id2 = 0;
843  std::vector<PointPair> inPoints;
844  for(i = 0; i < niso; ++i)
845  {
846  for(id1 = 0; id1 < iso[i]->m_nvert; ++id1)
847  {
848  inPoints.push_back(PointPair(BPoint( iso[i]->m_x[id1],
849  iso[i]->m_y[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);
853  id2++;
854  }
855  }
856 
857 
858  if(verbose)
859  {
860  cout << "Process building tree ..." << endl;
861  }
862 
863  //Build tree
864  bgi::rtree<PointPair, bgi::rstar<16> > rtree;
865  rtree.insert(inPoints.begin(), inPoints.end());
866 
867  //Find neipghbours
868  int unique_index = 0;
869  int prog=0;
870  for(i = 0; i < m_nvert; ++i)
871  {
872  if(verbose)
873  {
874  prog = LibUtilities::PrintProgressbar(i,m_nvert,
875  "Nearest verts",prog);
876  }
877 
878  BPoint queryPoint = inPoints[i].first;
879 
880 
881  // check to see if point has been already reset to lower than
882  // unique value
883  if(global_to_unique_map[i] < unique_index) // do nothing
884  {
885  }
886  else
887  {
888 
889  // find nearest 10 points within the distance box
890  std::vector<PointPair> result;
891  rtree.query(bgi::nearest(queryPoint, 10), std::back_inserter(result));
892 
893  //see if any values have unique value already
894  set<int> samept;
896  int new_index = -1;
897  for(id1 = 0; id1 < result.size(); ++id1)
898  {
899  NekDouble dist = bg::distance(queryPoint, result[id1].first);
900  if(dist*dist < NekConstants::kNekZeroTol) // same point
901  {
902  id2 = result[id1].second;
903  samept.insert(id2);
904 
905  if(global_to_unique_map[id2] <unique_index)
906  {
907  new_index = global_to_unique_map[id2];
908  }
909  }
910  }
911  if(new_index == -1)
912  {
913  new_index = unique_index;
914  unique_index++;
915  }
916 
917  // reset all same values to new_index
918  global_to_unique_map[i] = new_index;
919  for(it = samept.begin(); it != samept.end(); ++it)
920  {
921  global_to_unique_map[*it] = new_index;
922  }
923  }
924  }
925 
926  if(verbose)
927  {
928  cout << endl;
929  }
930 
931  m_nvert = unique_index;
932 
933  nelmt = 0;
934  // reset m_vid;
935  int cnt = 0;
936  for(n = 0; n < niso; ++n)
937  {
938  for(i = 0; i < iso[n]->m_ntris; ++i,nelmt++)
939  {
940  for(j=0; j < 3;++j)
941  {
942  m_vid[3*nelmt+j] = global_to_unique_map[iso[n]->m_vid[3*i+j]+cnt];
943  }
944  }
945  cnt += iso[n]->m_nvert;
946  }
947 
948  m_ntris = nelmt;
949 
950  m_x.resize(m_nvert);
951  m_y.resize(m_nvert);
952  m_z.resize(m_nvert);
953 
954  m_fields.resize(iso[0]->m_fields.size());
955  for(i = 0; i < iso[0]->m_fields.size(); ++i)
956  {
957  m_fields[i].resize(m_nvert);
958  }
959 
960  for(n = 0; n < global_to_unique_map.size(); ++n)
961  {
962 
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);
966 
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)
970  {
971  m_fields[j][global_to_unique_map[n]] = iso[isoid]->
972  m_fields[j][ptid];
973  }
974  }
975 }
976 
977 
978 // define == if point is within 1e-8
979 bool operator == (const IsoVertex& x, const IsoVertex& y)
980 {
981  return ((x.m_x-y.m_x)*(x.m_x-y.m_x) + (x.m_y-y.m_y)*(x.m_y-y.m_y) +
982  (x.m_z-y.m_z)*(x.m_z-y.m_z) < NekConstants::kNekZeroTol)? true:false;
983 }
984 
985 // define != if point is outside 1e-8
986 bool operator != (const IsoVertex& x, const IsoVertex& y)
987 {
988  return ((x.m_x-y.m_x)*(x.m_x-y.m_x) + (x.m_y-y.m_y)*(x.m_y-y.m_y) +
989  (x.m_z-y.m_z)*(x.m_z-y.m_z) < NekConstants::kNekZeroTol)? 0:1;
990 }
991 
992 void Iso::Smooth(int n_iter, NekDouble lambda, NekDouble mu)
993 {
994  int iter,i,j,k;
995  NekDouble del_v[3];
996  NekDouble w;
997  Array<OneD, NekDouble> xtemp, ytemp, ztemp;
998  vector< vector<int> > adj,vertcon;
999  vector< vector<NekDouble > > wght;
1002 
1003  // determine elements around each vertex
1004  vertcon.resize(m_nvert);
1005  for(i = 0; i < m_ntris; ++i)
1006  {
1007  for(j = 0; j < 3; ++j)
1008  {
1009  vertcon[m_vid[3*i+j]].push_back(i);
1010  }
1011  }
1012 
1013  // determine vertices and weights around each vertex
1014  adj.resize(m_nvert);
1015  wght.resize(m_nvert);
1016 
1017  for(i =0; i < m_nvert; ++i)
1018  {
1019  // loop over surrounding elements
1020  for(ipt = vertcon[i].begin(); ipt != vertcon[i].end(); ++ipt)
1021  {
1022  for(j = 0; j < 3; ++j)
1023  {
1024  // make sure not adding own vertex
1025  if(m_vid[3*(*ipt)+j] != i)
1026  {
1027  // check to see if vertex has already been added
1028  for(k = 0; k < adj[i].size(); ++k)
1029  {
1030  if(adj[i][k] == m_vid[3*(*ipt)+j])
1031  {
1032  break;
1033  }
1034  }
1035 
1036  if(k == adj[i].size())
1037  {
1038  adj[i].push_back(m_vid[3*(*ipt)+j]);
1039  }
1040  }
1041  }
1042  }
1043  }
1044 
1045  // Currently set weights up as even distribution
1046  for(i =0; i < m_nvert; ++i)
1047  {
1048  w = 1.0/((NekDouble)adj[i].size());
1049  for(j = 0; j < adj[i].size(); ++j)
1050  {
1051  wght[i].push_back(w);
1052  }
1053  }
1054 
1055 
1059 
1060  // smooth each point
1061  for (iter=0;iter<n_iter;iter++)
1062  {
1063  // compute first weighted average
1064  for(i=0;i< m_nvert;++i)
1065  {
1066  del_v[0] = del_v[1] = del_v[2] = 0.0;
1067  for(j = 0; j < adj[i].size(); ++j)
1068  {
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];
1072  }
1073 
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;
1077  }
1078 
1079  // compute second weighted average
1080  for(i=0;i< m_nvert;++i)
1081  {
1082  del_v[0] = del_v[1] = del_v[2] = 0;
1083  for(j = 0; j < adj[i].size(); ++j)
1084  {
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];
1088  }
1089 
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;
1093  }
1094  }
1095 }
1096 
1097 
1098 void Iso::SeparateRegions(vector<IsoSharedPtr> &sep_iso, int minsize, bool verbose)
1099 {
1100  int i,j,k,id;
1101  Array<OneD, vector<int> >vertcon(m_nvert);
1103  list<int> tocheck;
1104  list<int>::iterator cid;
1105 
1106  Array<OneD, bool> viddone(m_nvert,false);
1107 
1108  // make list of connecting tris around each vertex
1109  for(i = 0; i < m_ntris; ++i)
1110  {
1111  for(j = 0; j < 3; ++j)
1112  {
1113  vertcon[m_vid[3*i+j]].push_back(i);
1114  }
1115  }
1116 
1117  Array<OneD, int> vidregion(m_nvert,-1);
1118 
1119  int nregions = -1;
1120 
1121 
1122  // check all points are assigned to a region
1123  int progcnt = -1;
1124  for(k = 0; k < m_nvert; ++k)
1125  {
1126  if(verbose)
1127  {
1128  progcnt = LibUtilities::PrintProgressbar(k,m_nvert,"Separating regions",progcnt);
1129  }
1130 
1131  if(vidregion[k] == -1)
1132  {
1133  vidregion[k] = ++nregions;
1134 
1135  // find all elmts around this.. vertex that need to be checked
1136  for(ipt = vertcon[k].begin(); ipt != vertcon[k].end(); ++ipt)
1137  {
1138  for(i = 0; i < 3; ++i)
1139  {
1140  if(vidregion[id = m_vid[3*(ipt[0])+i]] == -1)
1141  {
1142  tocheck.push_back(id);
1143  vidregion[id] = nregions;
1144  }
1145  }
1146  }
1147  viddone[k] = 1;
1148 
1149  // check all other neighbouring vertices
1150  while(tocheck.size())
1151  {
1152  cid = tocheck.begin();
1153  while(cid != tocheck.end())
1154  {
1155  if(!viddone[*cid])
1156  {
1157  for(ipt = vertcon[*cid].begin(); ipt != vertcon[*cid].end(); ++ipt)
1158  {
1159  for(i = 0; i < 3; ++i)
1160  {
1161  if(vidregion[id = m_vid[3*(ipt[0])+i]] == -1)
1162  {
1163  tocheck.push_back(id);
1164  vidregion[id] = nregions;
1165  }
1166  }
1167  }
1168  viddone[*cid] = 1;
1169  ++cid;
1170  tocheck.pop_front();
1171  }
1172  }
1173  }
1174  }
1175  }
1176  nregions++;
1177 
1178 
1179  Array<OneD, int> nvert(nregions,0);
1180  Array<OneD, int> nelmt(nregions,0);
1181 
1182  // count nverts
1183  for(i = 0; i < m_nvert; ++i)
1184  {
1185  nvert[vidregion[i]] +=1;
1186  }
1187 
1188  // count nelmts
1189  for(i = 0; i < m_ntris; ++i)
1190  {
1191  nelmt[vidregion[m_vid[3*i]]] +=1;
1192  }
1193 
1194  Array<OneD, int> vidmap(m_nvert);
1195  // generate new list of isocontour
1196  for(int n = 0; n < nregions; ++n)
1197  {
1198  if(nelmt[n] > minsize)
1199  {
1200  int nfields = m_fields.size();
1202 
1203  iso->m_ntris = nelmt[n];
1204  iso->m_vid = Array<OneD, int>(3*nelmt[n]);
1205 
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]);
1210 
1211  iso->m_fields.resize(nfields);
1212  for(i = 0; i < nfields; ++i)
1213  {
1214  iso->m_fields[i].resize(nvert[n]);
1215  }
1216 
1217 
1218  int cnt = 0;
1219  // generate vid map;
1220  Vmath::Zero(m_nvert,vidmap,1);
1221  for(i = 0; i < m_nvert; ++i)
1222  {
1223  if(vidregion[i] == n)
1224  {
1225  vidmap[i] = cnt++;
1226  }
1227  }
1228 
1229  cnt = 0;
1230  for(i = 0; i < m_ntris; ++i)
1231  {
1232  if(vidregion[m_vid[3*i]] == n)
1233  {
1234  for(j = 0; j < 3; ++j)
1235  {
1236  iso->m_vid[3*cnt+j] = vidmap[m_vid[3*i+j]];
1237 
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]];
1241 
1242  for(k = 0; k < nfields; ++k)
1243  {
1244  iso->m_fields[k][vidmap[m_vid[3*i+j]]] = m_fields[k][m_vid[3*i+j]];
1245  }
1246  }
1247  cnt++;
1248  }
1249  }
1250 
1251  WARNINGL0(cnt == nelmt[n],"Number of elements do not match");
1252  sep_iso.push_back(iso);
1253  }
1254  }
1255 }
1256 
1257 }
1258 }
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
bool operator==(TriFaceIDs const &p1, TriFaceIDs const &p2)
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
Definition: Progressbar.hpp:69
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Represents a command-line configuration option.
STL namespace.
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)
void Stop()
Definition: Timer.cpp:62
static const NekDouble kNekZeroTol
boost::shared_ptr< Iso > IsoSharedPtr
This processing module interpolates one field to another.
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:767
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:204
static std::string npts
Definition: InputFld.cpp:43
void GlobalCondense(vector< boost::shared_ptr< Iso > > &iso, bool verbose)
double NekDouble
void Smooth(int n_iter, NekDouble lambda, NekDouble mu)
This class defines evaluator of analytic (symbolic) mathematical expressions. Expressions are allowed...
vector< NekDouble > m_z
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
vector< NekDouble > m_x
Array< OneD, int > m_vid
void SeparateRegions(vector< boost::shared_ptr< Iso > > &iso, int minsize, bool verbose)
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
void Start()
Definition: Timer.cpp:51
vector< IsoSharedPtr > ExtractContour(const int fieldid, const NekDouble val)
vector< NekDouble > m_y
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:108
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215