Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 using namespace std;
38 
39 #include "ProcessIsoContour.h"
40 
44 #include <boost/math/special_functions/fpclassify.hpp>
45 
46 namespace Nektar
47 {
48 namespace Utilities
49 {
50 
51 ModuleKey ProcessIsoContour::className =
53  ModuleKey(eProcessModule, "isocontour"),
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 "
58  "smoothing");
59 
60 ProcessIsoContour::ProcessIsoContour(FieldSharedPtr f) :
62 {
63 
64  m_config["fieldstr"] = ConfigOption(false, "NotSet",
65  "string of isocontour to be extracted");
66  m_config["fieldname"] = ConfigOption(false, "isocon",
67  "name for isocontour if fieldstr "
68  "specified, default is isocon");
69 
70  m_config["fieldid"] = ConfigOption(false, "NotSet",
71  "field id to extract");
72 
73  m_config["fieldvalue"] = ConfigOption(false, "NotSet",
74  "field value to extract");
75 
76  m_config["globalcondense"] = ConfigOption(true, "NotSet",
77  "Globally condense contour to unique "
78  "values");
79 
80  m_config["smooth"] = ConfigOption(true, "NotSet",
81  "Smooth isocontour (might require "
82  "globalcondense)");
83 
84  m_config["smoothiter"] = ConfigOption(false, "100",
85  "Number of smoothing cycle, default = "
86  "100");
87 
88  m_config["smoothposdiffusion"] = ConfigOption(false,"0.5",
89  "Postive diffusion coefficient "
90  "(0 < lambda < 1), default = 0.5");
91 
92  m_config["smoothnegdiffusion"] = ConfigOption(false,"0.495",
93  "Negative diffusion coefficient "
94  "(0 < mu < 1), default = 0.495");
95 
96  m_config["removesmallcontour"] = ConfigOption(false,"0",
97  "Remove contours with less than specified number of triangles."
98  "Only valid with GlobalCondense or Smooth options.");
99 
100 
101 }
102 
104 {
105 }
106 
107 void ProcessIsoContour::Process(po::variables_map &vm)
108 {
109  Timer timer;
110  int rank = m_f->m_comm->GetRank();
111 
112  if(m_f->m_verbose)
113  {
114  if(rank == 0)
115  {
116  cout << "Process Contour extraction..." << endl;
117  timer.Start();
118  }
119  }
120 
121  vector<IsoSharedPtr> iso;
122 
123  if(m_f->m_fieldPts.get()) // assume we have read .dat file to directly input dat file.
124  {
126  }
127  else // extract isocontour from field
128  {
129  if(m_f->m_exp.size() == 0)
130  {
131  return;
132  }
133 
134  // extract all fields to equi-spaced
136 
137  int fieldid;
138  NekDouble value;
139 
140  if(m_config["fieldstr"].m_beenSet) //generate field of interest
141  {
142  fieldid = m_f->m_fieldPts->GetNFields();
143 
144  Array<OneD, NekDouble> pts(m_f->m_fieldPts->GetNpoints());
145 
146  // evaluate new function
148  string varstr = "x y z";
149  vector<Array<OneD, const NekDouble> > interpfields;
150 
151  for(int i = 0; i < m_f->m_fieldPts->GetDim(); ++i)
152  {
153  interpfields.push_back(m_f->m_fieldPts->GetPts(i));
154  }
155  for(int i = 0; i < m_f->m_fieldPts->GetNFields(); ++i)
156  {
157  varstr += " " + m_f->m_fieldPts->GetFieldName(i);
158  interpfields.push_back(m_f->m_fieldPts->GetPts(i+3));
159  }
160 
161  int ExprId = -1;
162  std::string str = m_config["fieldstr"].as<string>();
163  ExprId = strEval.DefineFunction(varstr.c_str(), str);
164 
165  strEval.Evaluate(ExprId, interpfields, pts);
166 
167  // set up field name if provided otherwise called "isocon" from default
168  string fieldName = m_config["fieldname"].as<string>();
169 
170  m_f->m_fieldPts->AddField(pts, fieldName);
171  }
172  else
173  {
174  ASSERTL0(m_config["fieldid"].as<string>() != "NotSet", "fieldid must be specified");
175  fieldid = m_config["fieldid"].as<int>();
176  }
177 
178  ASSERTL0(m_config["fieldvalue"].as<string>() != "NotSet", "fieldvalue must be specified");
179  value = m_config["fieldvalue"].as<NekDouble>();
180 
181  iso = ExtractContour(fieldid,value);
182 
183  }
184 
185  // Process isocontour
186  bool smoothing = m_config["smooth"].m_beenSet;
187  bool globalcondense = m_config["globalcondense"].m_beenSet;
188  if(globalcondense)
189  {
190  int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
192 
193  g_iso->globalcondense(iso,m_f->m_verbose);
194 
195 
196  iso.clear();
197  iso.push_back(g_iso);
198  }
199 
200  if(smoothing)
201  {
202  int niter = m_config["smoothiter"].as<int>();
203  NekDouble lambda = m_config["smoothposdiffusion"].as<NekDouble>();
204  NekDouble mu = m_config["smoothnegdiffusion"].as<NekDouble>();
205  for(int i =0 ; i < iso.size(); ++i)
206  {
207  iso[i]->smooth(niter,lambda,-mu);
208  }
209  }
210 
211  int mincontour = 0;
212  if((mincontour = m_config["removesmallcontour"].as<int>()))
213  {
214  vector<IsoSharedPtr> new_iso;
215 
216  if(rank == 0)
217  {
218  cout << "Identifying separate regions [." << flush ;
219  }
220  for(int i =0 ; i < iso.size(); ++i)
221  {
222  iso[i]->separate_regions(new_iso,mincontour,m_f->m_verbose);
223  }
224 
225  if(rank == 0)
226  {
227  cout << "]" << endl << flush ;
228  }
229 
230  // reset iso to new_iso;
231  iso = new_iso;
232  }
233 
234  ResetFieldPts(iso);
235 
236 
237  if(m_f->m_verbose)
238  {
239  if(rank == 0)
240  {
241  timer.Stop();
242  NekDouble cpuTime = timer.TimePerTest(1);
243 
244  stringstream ss;
245  ss << cpuTime << "s";
246  cout << "Process Isocontour CPU Time: " << setw(8) << left
247  << ss.str() << endl;
248  cpuTime = 0.0;
249  }
250  }
251 }
252 
253 
254 /*********************/
255 /* Function TwoPairs */
256 /*********************/
257 
261  int &pr)
262 /* The logic of the following SWITCH statement may only be
263  understood with a "peusdo truth-table" */
264 {
265  if (((cx[0]-cx[1])==0.0)&&
266  ((cy[0]-cy[1])==0.0)&&
267  ((cz[0]-cz[1])==0.0))
268  {
269  if (((cx[2]-cx[3])==0.0)&&
270  ((cy[2]-cy[3])==0.0)&&
271  ((cz[2]-cz[3])==0.0))
272  {
273  pr=4;
274  }
275  else
276  {
277  pr=3;
278  }
279  }
280  else
281  {
282  pr=1;
283  }
284 }
285 
286 
287 /*************************/
288 /* Function ThreeSimilar */
289 /*************************/
290 void ThreeSimilar (const int i, const int j, const int k,
291  int &pr, int &ps)
292 /* The logic of the following SWITCH statement may be
293  understood with a "peusdo truth-table" */
294 {
295  switch (i + j + k)
296  {
297  case (3):
298  pr = 3;
299  ps = 4;
300  break;
301  case (4):
302  pr = 2;
303  ps = 4;
304  break;
305  case (5):
306  if (j == 1)
307  {
308  pr = 2;
309  ps = 3;
310  }
311  else
312  {
313  pr = 1;
314  ps = 4;
315  }
316  break;
317  case (6):
318  if (i == 0)
319  {
320  pr = 1;
321  ps = 3;
322  }
323  else
324  {
325  pr = 0;
326  ps = 4;
327  }
328  break;
329  case (7):
330  if (i == 0)
331  {
332  pr = 1;
333  ps = 2;
334  }
335  else
336  {
337  pr = 0;
338  ps = 3;
339  }
340  break;
341  case (8):
342  pr = 0;
343  ps = 2;
344  break;
345  case (9):
346  pr = 0;
347  ps = 1;
348  break;
349  default:
350  ASSERTL0(false,"Error in 5-point triangulation in ThreeSimilar");
351  break;
352  }
353 }
354 
356  const int fieldid,
357  const NekDouble val)
358 {
359  vector<IsoSharedPtr> returnval;
360 
361  int coordim = m_f->m_fieldPts->GetDim();
362  int nfields = m_f->m_fieldPts->GetNFields() + coordim;
363 
364  ASSERTL0(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock,
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);
370 
371  Array<OneD, NekDouble> c = fields[coordim + fieldid];
372 
373  int i, j, k, ii, jj, kk, r, s, n, counter, boolean;
374  Array<OneD, Array<OneD, NekDouble> > intfields(nfields);
375  intfields[0] = Array<OneD, NekDouble>(5*nfields);
376  for(i = 1; i < nfields; ++i)
377  {
378  intfields[i] = intfields[i-1] + 5;
379  }
380  Array<OneD, NekDouble> cx = intfields[0];
381  Array<OneD, NekDouble> cy = intfields[1];
382  Array<OneD, NekDouble> cz = intfields[2];
383 
384  vector<Array<OneD, int> > ptsConn;
385  m_f->m_fieldPts->GetConnectivity(ptsConn);
386  for(int zone = 0; zone < ptsConn.size(); ++zone)
387  {
388  IsoSharedPtr iso;
389 
390  iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-3);
391 
392  int nelmt = ptsConn[zone].num_elements()
393  /(coordim+1);
394 
395  Array<OneD, int> conn = ptsConn[zone];
396 
397  for (n = 0, i = 0; i < nelmt; ++i)
398  {
399  // check to see if val is between vertex values
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))))
404  {
405 
406  // loop over all edges and interpolate if
407  // contour is between vertex values
408  for (counter = 0, j=0; j<=2; j++)
409  {
410  for (k=j+1; k<=3; k++)
411  {
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]])))
416  {
417  // linear interpolation of fields
418  // (and coords).
419  NekDouble cj = c[conn[i*4+j]];
420  NekDouble ck = c[conn[i*4+k]];
421  NekDouble factor = (val-cj)/(ck-cj);
422 
423  if(fabs(cj-ck) > 1e-12)
424  {
425  // interpolate coordinates and fields
426  for(int f = 0; f < nfields; ++f)
427  {
428  if(counter == 5)
429  {
430  ASSERTL0(false,"Counter is 5");
431  }
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]]);
436  }
437  ++counter;
438  }
439  }
440  }
441  }
442 
443  switch(counter)
444  {
445  case 3:
446  n+=1;
447  iso->resize_fields(3*n);
448 
449  for(j = 0; j < 3; ++j)
450  {
451  iso->set_fields(3*(n-1)+j,intfields,j);
452  }
453  break;
454  case 4:
455  n+=2;
456  iso->resize_fields(3*n);
457 
458  for(j = 0; j < 3; ++j)
459  {
460  iso->set_fields(3*(n-2)+j,intfields,j);
461  iso->set_fields(3*(n-1)+j,intfields,j+1);
462  }
463  break;
464  case 5:
465  n+=1;
466  iso->resize_fields(3*n);
467 
468  boolean=0;
469  for (ii=0;ii<=2;ii++)
470  {
471  for (jj=ii+1;jj<=3;jj++)
472  {
473  for (kk=jj+1;kk<=4;kk++)
474  {
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)))
481  {
482  boolean+=1;
483  ThreeSimilar (ii,jj,kk,r,s);
484 
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);
488  }
489  else
490  {
491  boolean+=0;
492  }
493  }
494  }
495  }
496 
497  if (boolean==0)
498  {
499  TwoPairs (cx,cy,cz,r);
500 
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);
504  }
505  break;
506  }
507  }
508  }
509  iso->set_ntris(n);
510 
511  // condense the information in this elemental extraction.
512  iso->condense();
513  returnval.push_back(iso);
514  }
515 
516  return returnval;
517 }
518 
519 // reset m_fieldPts with values from iso;
520 void ProcessIsoContour::ResetFieldPts(vector<IsoSharedPtr> &iso)
521 {
522  int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
523 
524  // set output to triangle block.
525  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
526 
527  Array<OneD, Array<OneD, NekDouble> > newfields(nfields);
528 
529  // count up number of points
530  int npts = 0;
531  for(int i =0; i < iso.size(); ++i)
532  {
533  npts += iso[i]->get_nvert();
534  }
535 
536  // set up coordinate in new field
537  newfields[0] = Array<OneD, NekDouble>(npts);
538  newfields[1] = Array<OneD, NekDouble>(npts);
539  newfields[2] = Array<OneD, NekDouble>(npts);
540 
541  int cnt = 0;
542  for(int i =0; i < iso.size(); ++i)
543  {
544  for(int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
545  {
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);
549  }
550  }
551 
552 
553  // set up fields
554  for(int f = 0; f < nfields-3; ++f)
555  {
556  newfields[f+3] = Array<OneD, NekDouble>(npts);
557 
558  cnt = 0;
559  for(int i =0; i < iso.size(); ++i)
560  {
561  for(int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
562  {
563  newfields[f+3][cnt] = iso[i]->get_fields(f,j);
564  }
565  }
566  }
567 
568  m_f->m_fieldPts->SetPts(newfields);
569 
570  // set up connectivity data.
571  vector<Array<OneD, int> > ptsConn;
572  m_f->m_fieldPts->GetConnectivity(ptsConn);
573  cnt = 0;
574  ptsConn.clear();
575  for(int i =0; i < iso.size(); ++i)
576  {
577  int ntris = iso[i]->get_ntris();
578  Array<OneD, int> conn(ntris*3);
579 
580  for(int j = 0; j < 3*ntris; ++j)
581  {
582  conn[j] = cnt + iso[i]->get_vid(j);
583  }
584  ptsConn.push_back(conn);
585  cnt += iso[i]->get_nvert();
586  }
587  m_f->m_fieldPts->SetConnectivity(ptsConn);
588 }
589 
590 // reset m_fieldPts with values from iso;
591 void ProcessIsoContour::SetupIsoFromFieldPts(vector<IsoSharedPtr> &isovec)
592 {
593  ASSERTL0(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock,
594  "Assume input is from ePtsTriBlock");
595 
596  // get information from PtsField
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);
603 
604 
605  int cnt = 0;
606  for(int c = 0; c < ptsConn.size(); ++c)
607  {
608  // set up single iso with all the information from PtsField
610 
611  int nelmt = 0;
612  nelmt = ptsConn[c].num_elements()/3;
613 
614  iso->set_ntris(nelmt);
615  iso->resize_vid(3*nelmt);
616 
617  // fill in connectivity values.
618  int nvert = 0;
619  for(int i = 0; i < ptsConn[c].num_elements(); ++i)
620  {
621  int cid = ptsConn[c][i]-cnt;
622  iso->set_vid(i,cid);
623  nvert = max(cid,nvert);
624  }
625  nvert++;
626 
627  iso->set_nvert(nvert);
628  iso->resize_fields(nvert);
629 
630  // fill in points values (including coordinates)
631  for(int i = 0; i < nvert; ++i)
632  {
633  iso->set_fields(i,fieldpts,i+cnt);
634  }
635  cnt += nvert;
636  isovec.push_back(iso);
637  }
638 
639 }
640 
641 
642 void Iso::condense(void)
643 {
644  register int i,j,cnt;
645  IsoVertex v;
646  vector<IsoVertex> vert;
648 
649  if(!m_ntris) return;
650 
651  if(m_condensed) return;
652  m_condensed = true;
653 
654  vert.reserve(m_ntris/6);
655 
657 
658  // fill first 3 points and initialise fields
659  v.m_fields.resize(m_fields.size());
660  for(cnt =0, i = 0; i < 3; ++i)
661  {
662  v.m_x = m_x[i];
663  v.m_y = m_y[i];
664  v.m_z = m_z[i];
665  for(int f = 0; f < m_fields.size(); ++f)
666  {
667  v.m_fields[f] = m_fields[f][i];
668  }
669  v.m_id = cnt;
670  vert.push_back(v);
671  m_vid[i] = v.m_id;
672  ++cnt;
673  }
674 
675  for(i = 1; i < m_ntris; ++i)
676  {
677  for(j = 0; j < 3; ++j)
678  {
679  v.m_x = m_x[3*i+j];
680  v.m_y = m_y[3*i+j];
681  v.m_z = m_z[3*i+j];
682 
683  pt = find(vert.begin(),vert.end(),v);
684  if(pt != vert.end())
685  {
686  m_vid[3*i+j] = pt[0].m_id;
687  }
688  else
689  {
690  v.m_id = cnt;
691 
692  for(int f = 0; f < m_fields.size(); ++f)
693  {
694  v.m_fields[f] = m_fields[f][3*i+j];
695  }
696 
697  vert.push_back(v);
698 
699  m_vid[3*i+j] = v.m_id;
700  ++cnt;
701  }
702  }
703  }
704 
705  // remove elements with multiple vertices
706  for(i = 0,cnt=0; i < m_ntris;)
707  {
708  if((m_vid[3*i] ==m_vid[3*i+1])||
709  (m_vid[3*i] ==m_vid[3*i+2])||
710  (m_vid[3*i+1]==m_vid[3*i+2]))
711  {
712  cnt++;
713  for(j = 3*i; j < 3*(m_ntris-1); ++j)
714  {
715  m_vid[j] = m_vid[j+3];
716  }
717  m_ntris--;
718  }
719  else
720  {
721  ++i;
722  }
723  }
724 
725  m_nvert = vert.size();
726 
727  m_x.resize(m_nvert);
728  m_y.resize(m_nvert);
729  m_z.resize(m_nvert);
730 
731  for(int f = 0; f < m_fields.size(); ++f)
732  {
733  m_fields[f].resize(m_nvert);
734  }
735 
736  for(i = 0; i < m_nvert; ++i)
737  {
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)
742  {
743  m_fields[f][i] = vert[i].m_fields[f];
744  }
745  }
746 }
747 
748 
750 
751 // define == if point is within 1e-4
752 bool operator == (const IsoVertex& x, const IsoVertex& y)
753 {
754  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) +
755  (x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? true:false;
756 }
757 
758 // define != if point is outside 1e-4
759 bool operator != (const IsoVertex& x, const IsoVertex& y)
760 {
761  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) +
762  (x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? 0:1;
763 }
764 
765 
767  NekDouble x2, NekDouble y2, NekDouble z2)
768 {
769  if((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2) < SQ_PNT_TOL)
770  {
771  return true;
772  }
773 
774  return false;
775 }
776 
777 void Iso::globalcondense(vector<IsoSharedPtr> &iso, bool verbose)
778 {
779  int i,j,n;
780  int nvert,nelmt;
781  int niso=iso.size();
782  int id1,id2;
784 
785  if(m_condensed) return;
786  m_condensed = true;
787 
788  vidmap = Array<OneD, Array<OneD, int> > (niso);
789 
790  m_ntris = 0;
791  for(i = 0; i < niso; ++i)
792  {
793  if(iso[i]->m_ntris)
794  {
795  m_ntris += iso[i]->m_ntris;
796  }
797  }
798 
800 
801  m_nvert = 0;
802  for(i = 0; i < niso; ++i)
803  {
804  if(iso[i]->m_ntris)
805  {
806  m_nvert += iso[i]->m_nvert;
807  }
808  }
809 
810  vector< vector<int> > isocon;
811  isocon.resize(niso);
812 
813  // identify which iso are connected by at least one point;
814  // find min x,y,z and max x,y,z and see if overlap to select
815  // which zones should be connected
816  {
817  vector<Array<OneD, NekDouble> > sph(niso);
818  Array<OneD, NekDouble> rng(6);
819  for(i = 0; i < niso; ++i)
820  {
821  sph[i] = Array<OneD, NekDouble>(4);
822 
823  // find max and min of isocontour
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];
827 
828  for(id1 = 1; id1 < iso[i]->m_nvert;++id1)
829  {
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]);
833 
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]);
837  }
838 
839  // centroid
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;
843 
844  // radius;
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]));
848  }
849 
850  for(i = 0; i < niso; ++i)
851  {
852  for(j = i; j < niso; ++j)
853  {
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]));
857 
858  // if centroid is closer than added radii
859  if(diff < sph[i][3] + sph[j][3])
860  {
861  isocon[i].push_back(j);
862  }
863  }
864  }
865 
866  }
867 
868 
869  for(i = 0; i < niso; ++i)
870  {
871  vidmap[i] = Array<OneD, int>(iso[i]->m_nvert,-1);
872  }
873  nvert = 0;
874  int cnt = 0;
875  // count up amount of checking to be done
876  NekDouble totiso = 0;
877  for(i = 0; i < niso; ++i)
878  {
879  totiso += isocon[i].size();
880  }
881 
882 
883  if(verbose)
884  {
885  cout << "Progress Bar totiso: " << totiso << endl;
886  }
887  for(i = 0; i < niso; ++i)
888  {
889  for(n = 0; n < isocon[i].size(); ++n, ++cnt)
890  {
891 
892  if(verbose && totiso >= 40)
893  {
894  LibUtilities::PrintProgressbar(cnt,totiso,"Condensing verts");
895  }
896 
897  int con = isocon[i][n];
898  for(id1 = 0; id1 < iso[i]->m_nvert; ++id1)
899  {
900 
901  if(verbose && totiso < 40)
902  {
903  LibUtilities::PrintProgressbar(id1,iso[i]->m_nvert,"isocon");
904  }
905 
906  int start = 0;
907  if(con == i)
908  {
909  start = id1+1;
910  }
911  for(id2 = start; id2 < iso[con]->m_nvert; ++id2)
912  {
913 
914  if((vidmap[con][id2] == -1)||(vidmap[i][id1] == -1))
915  {
916  if(same(iso[i]->m_x[id1], iso[i]->m_y[id1],
917  iso[i]->m_z[id1], iso[con]->m_x[id2],
918  iso[con]->m_y[id2],iso[con]->m_z[id2]))
919  {
920  if((vidmap[i][id1] == -1) &&
921  (vidmap[con][id2] != -1))
922  {
923  vidmap[i][id1] = vidmap[con][id2];
924  }
925  else if((vidmap[con][id2] == -1) &&
926  (vidmap[i][id1] != -1))
927  {
928  vidmap[con][id2] = vidmap[i][id1];
929  }
930  else if((vidmap[con][id2] == -1) &&
931  (vidmap[i][id1] == -1))
932  {
933  vidmap[i][id1] = vidmap[con][id2] = nvert++;
934  }
935  }
936  }
937  }
938  }
939  }
940 
941  for(id1 = 0; id1 < iso[i]->m_nvert;++id1)
942  {
943  if(vidmap[i][id1] == -1)
944  {
945  vidmap[i][id1] = nvert++;
946  }
947  }
948  }
949  m_nvert = nvert;
950 
951  nelmt = 0;
952  // reset m_vid;
953  for(n = 0; n < niso; ++n)
954  {
955  for(i = 0; i < iso[n]->m_ntris; ++i,nelmt++)
956  {
957  for(j=0; j < 3;++j)
958  {
959  m_vid[3*nelmt+j] = vidmap[n][iso[n]->m_vid[3*i+j]];
960  }
961  }
962  }
963 
964  m_ntris = nelmt;
965 
966  m_x.resize(m_nvert);
967  m_y.resize(m_nvert);
968  m_z.resize(m_nvert);
969 
970  m_fields.resize(iso[0]->m_fields.size());
971  for(i = 0; i < iso[0]->m_fields.size(); ++i)
972  {
973  m_fields[i].resize(m_nvert);
974  }
975 
976  // reset coordinate and fields.
977  for(n = 0; n < niso; ++n)
978  {
979  for(i = 0; i < iso[n]->m_nvert; ++i)
980  {
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];
984 
985  for(j = 0; j < m_fields.size(); ++j)
986  {
987  m_fields[j][vidmap[n][i]] = iso[n]->m_fields[j][i];
988  }
989  }
990  }
991  cout << endl;
992 }
993 
994 void Iso::smooth(int n_iter, NekDouble lambda, NekDouble mu)
995 {
996  int iter,i,j;
997  NekDouble del_v[3];
998  NekDouble w;
999  Array<OneD, NekDouble> xtemp, ytemp, ztemp;
1000  vector< vector<int> > adj,vertcon;
1003 
1004  // determine elements around each vertex
1005  vertcon.resize(m_nvert);
1006  for(i = 0; i < m_ntris; ++i)
1007  {
1008  for(j = 0; j < 3; ++j)
1009  {
1010  vertcon[m_vid[3*i+j]].push_back(i);
1011  }
1012  }
1013 
1014  // determine vertices around each vertex
1015  adj.resize(m_nvert);
1016 
1017  for(i =0; i < m_nvert; ++i)
1018  {
1019  for(ipt = vertcon[i].begin(); ipt != vertcon[i].end(); ++ipt)
1020  {
1021  for(j = 0; j < 3; ++j)
1022  {
1023  // make sure not adding own vertex
1024  if(m_vid[3*(*ipt)+j] != i)
1025  {
1026  // check to see if vertex has already been added
1027  for(iad = adj[i].begin(); iad != adj[i].end();++iad)
1028  {
1029  if(*iad == (m_vid[3*(*ipt)+j])) break;
1030  }
1031 
1032  if(iad == adj[i].end())
1033  {
1034  adj[i].push_back(m_vid[3*(*ipt)+j]);
1035  }
1036  }
1037  }
1038  }
1039  }
1040 
1044 
1045  // smooth each point
1046  for (iter=0;iter<n_iter;iter++)
1047  {
1048  // compute first weighted average
1049  for(i=0;i< m_nvert;++i)
1050  {
1051  w = 1.0/(NekDouble)(adj[i].size());
1052 
1053  del_v[0] = del_v[1] = del_v[2] = 0.0;
1054 
1055  for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
1056  {
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;
1060  }
1061 
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;
1065  }
1066 
1067  // compute second weighted average
1068  for(i=0;i< m_nvert;++i)
1069  {
1070 
1071  w = 1.0/(NekDouble)(adj[i].size());
1072  del_v[0] = del_v[1] = del_v[2] = 0;
1073 
1074  for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
1075  {
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;
1079  }
1080 
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;
1084  }
1085  }
1086 }
1087 
1088 
1089  void Iso::separate_regions(vector<IsoSharedPtr> &sep_iso, int minsize, bool verbose)
1090  {
1091  int i,j,k,id;
1092  Array<OneD, vector<int> >vertcon(m_nvert);
1094  list<int> tocheck;
1095  list<int>::iterator cid;
1096 
1097  Array<OneD, bool> viddone(m_nvert,false);
1098 
1099  // make list of connecting tris around each vertex
1100  for(i = 0; i < m_ntris; ++i)
1101  {
1102  for(j = 0; j < 3; ++j)
1103  {
1104  vertcon[m_vid[3*i+j]].push_back(i);
1105  }
1106  }
1107 
1108  Array<OneD, int> vidregion(m_nvert,-1);
1109 
1110  int nregions = -1;
1111 
1112 
1113  // check all points are assigned to a region
1114  for(k = 0; k < m_nvert; ++k)
1115  {
1116  if(verbose)
1117  {
1118  LibUtilities::PrintProgressbar(k,m_nvert,"Separating regions");
1119  }
1120 
1121  if(vidregion[k] == -1)
1122  {
1123  vidregion[k] = ++nregions;
1124 
1125  // find all elmts around this.. vertex that need to be checked
1126  for(ipt = vertcon[k].begin(); ipt != vertcon[k].end(); ++ipt)
1127  {
1128  for(i = 0; i < 3; ++i)
1129  {
1130  if(vidregion[id = m_vid[3*(ipt[0])+i]] == -1)
1131  {
1132  tocheck.push_back(id);
1133  vidregion[id] = nregions;
1134  }
1135  }
1136  }
1137  viddone[k] = 1;
1138 
1139  // check all other neighbouring vertices
1140  while(tocheck.size())
1141  {
1142  cid = tocheck.begin();
1143  while(cid != tocheck.end())
1144  {
1145  if(!viddone[*cid])
1146  {
1147  for(ipt = vertcon[*cid].begin(); ipt != vertcon[*cid].end(); ++ipt)
1148  {
1149  for(i = 0; i < 3; ++i)
1150  {
1151  if(vidregion[id = m_vid[3*(ipt[0])+i]] == -1)
1152  {
1153  tocheck.push_back(id);
1154  vidregion[id] = nregions;
1155  }
1156  }
1157  }
1158  viddone[*cid] = 1;
1159  ++cid;
1160  tocheck.pop_front();
1161  }
1162  }
1163  }
1164  }
1165  }
1166  nregions++;
1167 
1168 
1169  Array<OneD, int> nvert(nregions,0);
1170  Array<OneD, int> nelmt(nregions,0);
1171 
1172  // count nverts
1173  for(i = 0; i < m_nvert; ++i)
1174  {
1175  nvert[vidregion[i]] +=1;
1176  }
1177 
1178  // count nelmts
1179  for(i = 0; i < m_ntris; ++i)
1180  {
1181  nelmt[vidregion[m_vid[3*i]]] +=1;
1182  }
1183 
1184  Array<OneD, int> vidmap(m_nvert);
1185  // generate new list of isocontour
1186  for(int n = 0; n < nregions; ++n)
1187  {
1188  if(nelmt[n] > minsize)
1189  {
1190  int nfields = m_fields.size();
1192 
1193  iso->m_ntris = nelmt[n];
1194  iso->m_vid = Array<OneD, int>(3*nelmt[n]);
1195 
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]);
1200 
1201  iso->m_fields.resize(nfields);
1202  for(i = 0; i < nfields; ++i)
1203  {
1204  iso->m_fields[i].resize(nvert[n]);
1205  }
1206 
1207 
1208  int cnt = 0;
1209  // generate vid map;
1210  Vmath::Zero(m_nvert,vidmap,1);
1211  for(i = 0; i < m_nvert; ++i)
1212  {
1213  if(vidregion[i] == n)
1214  {
1215  vidmap[i] = cnt++;
1216  }
1217  }
1218 
1219  cnt = 0;
1220  for(i = 0; i < m_ntris; ++i)
1221  {
1222  if(vidregion[m_vid[3*i]] == n)
1223  {
1224  for(j = 0; j < 3; ++j)
1225  {
1226  iso->m_vid[3*cnt+j] = vidmap[m_vid[3*i+j]];
1227 
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]];
1231 
1232  for(k = 0; k < nfields; ++k)
1233  {
1234  iso->m_fields[k][vidmap[m_vid[3*i+j]]] = m_fields[k][m_vid[3*i+j]];
1235  }
1236  }
1237  cnt++;
1238  }
1239  }
1240 
1241  WARNINGL0(cnt == nelmt[n],"Number of elements do not match");
1242  sep_iso.push_back(iso);
1243  }
1244  }
1245  }
1246 }
1247 }
1248 
1249 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void Process()=0
vector< NekDouble > m_x
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
FieldSharedPtr m_f
Field object.
void separate_regions(vector< boost::shared_ptr< Iso > > &iso, int minsize, bool verbose)
void Stop()
Definition: Timer.cpp:62
void smooth(int n_iter, NekDouble lambda, NekDouble mu)
boost::shared_ptr< Iso > IsoSharedPtr
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:167
bool operator!=(const IsoVertex &x, const IsoVertex &y)
static std::string npts
Definition: InputFld.cpp:43
void SetupIsoFromFieldPts(vector< IsoSharedPtr > &isovec)
void TwoPairs(Array< OneD, NekDouble > &cx, Array< OneD, NekDouble > &cy, Array< OneD, NekDouble > &cz, int &pr)
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
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< NekDouble > m_z
vector< IsoSharedPtr > ExtractContour(const int fieldid, const NekDouble val)
Array< OneD, int > m_vid
vector< NekDouble > m_y
void globalcondense(vector< boost::shared_ptr< Iso > > &iso, bool verbose)
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
void ResetFieldPts(vector< IsoSharedPtr > &iso)
void Start()
Definition: Timer.cpp:51
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.
Definition: Vmath.cpp:359
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:191
This processing module interpolates one field to another.
void PrintProgressbar(const int position, const int goal, const string message)
Prints a progressbar.
Definition: Progressbar.hpp:69
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215