Nektar++
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 
43 #include <boost/math/special_functions/fpclassify.hpp>
44 
45 namespace Nektar
46 {
47 namespace Utilities
48 {
49 
50 ModuleKey ProcessIsoContour::className =
52  ModuleKey(eProcessModule, "isocontour"),
53  ProcessIsoContour::create,
54  "Extract an isocontour of fieldid variable and at "
55  "value fieldvalue, Optionally fieldstr can be "
56  "specified for a string defiition or smooth for "
57  "smoothing");
58 
59 ProcessIsoContour::ProcessIsoContour(FieldSharedPtr f) :
61 {
62 
63  m_config["fieldstr"] = ConfigOption(false, "NotSet",
64  "string of isocontour to be extracted");
65  m_config["fieldname"] = ConfigOption(false, "isocon",
66  "name for isocontour if fieldstr "
67  "specified, default is isocon");
68 
69  m_config["fieldid"] = ConfigOption(false, "NotSet",
70  "field id to extract");
71  m_config["fieldvalue"] = ConfigOption(false, "NotSet",
72  "field value to extract");
73 
74  m_config["globalcondense"] = ConfigOption(true, "NotSet",
75  "Globally condense contour to unique "
76  "values");
77 
78  m_config["smooth"] = ConfigOption(true, "NotSet",
79  "Smooth isocontour (implies global "
80  "condense)");
81  m_config["smoothiter"] = ConfigOption(false, "100",
82  "Number of smoothing cycle, default = "
83  "100");
84  m_config["smoothposdiffusion"] = ConfigOption(false,"0.5",
85  "Postive diffusion coefficient "
86  "(0 < lambda < 1), default = 0.5");
87  m_config["smoothnegdiffusion"] = ConfigOption(false,"0.495",
88  "Negative diffusion coefficient "
89  "(0 < mu < 1), default = 0.495");
90 
91  m_config["removesmallcontour"] = ConfigOption(false,"0",
92  "Remove contours with less than specified number of triangles."
93  "Only valid with GlobalCondense or Smooth options.");
94 
95 
96 }
97 
99 {
100 }
101 
102 void ProcessIsoContour::Process(po::variables_map &vm)
103 {
104  vector<IsoSharedPtr> iso;
105 
106  // extract all fields to equi-spaced
108 
109  int fieldid;
110  NekDouble value;
111 
112  if(m_config["fieldstr"].m_beenSet) //generate field of interest
113  {
114  fieldid = m_f->m_fieldPts->GetNFields();
115 
116  Array<OneD, NekDouble> pts(m_f->m_fieldPts->GetNpoints());
117 
118  // evaluate new function
120  string varstr = "x y z";
121  vector<Array<OneD, const NekDouble> > interpfields;
122 
123  for(int i = 0; i < m_f->m_fieldPts->GetDim(); ++i)
124  {
125  interpfields.push_back(m_f->m_fieldPts->GetPts(i));
126  }
127  for(int i = 0; i < m_f->m_fieldPts->GetNFields(); ++i)
128  {
129  varstr += " " + m_f->m_fieldPts->GetFieldName(i);
130  interpfields.push_back(m_f->m_fieldPts->GetPts(i+3));
131  }
132 
133  int ExprId = -1;
134  std::string str = m_config["fieldstr"].as<string>();
135  ExprId = strEval.DefineFunction(varstr.c_str(), str);
136 
137  strEval.Evaluate(ExprId, interpfields, pts);
138 
139  // set up field name if provided otherwise called "isocon" from default
140  string fieldName = m_config["fieldname"].as<string>();
141 
142  m_f->m_fieldPts->AddField(pts, fieldName);
143  }
144  else
145  {
146  ASSERTL0(m_config["fieldid"].as<string>() != "NotSet", "fieldid must be specified");
147  fieldid = m_config["fieldid"].as<int>();
148  }
149 
150  ASSERTL0(m_config["fieldvalue"].as<string>() != "NotSet", "fieldvalue must be specified");
151  value = m_config["fieldvalue"].as<NekDouble>();
152 
153  iso = ExtractContour(fieldid,value);
154  bool smoothing = m_config["smooth"].m_beenSet;
155  bool globalcondense = m_config["globalcondense"].m_beenSet;
156  if(smoothing||globalcondense)
157  {
158  vector<IsoSharedPtr> glob_iso;
159  int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
161  int mincontour = 0;
162 
163  g_iso->globalcondense(iso);
164 
165  if(smoothing)
166  {
167  int niter = m_config["smoothiter"].as<int>();
168  NekDouble lambda = m_config["smoothposdiffusion"].as<NekDouble>();
169  NekDouble mu = m_config["smoothnegdiffusion"].as<NekDouble>();
170  g_iso->smooth(niter,lambda,-mu);
171  }
172 
173 
174  if((mincontour = m_config["removesmallcontour"].as<int>()))
175  {
176  g_iso->separate_regions(glob_iso,mincontour);
177  }
178  else
179  {
180  glob_iso.push_back(g_iso);
181  }
182  ResetFieldPts(glob_iso);
183  }
184  else
185  {
186  ResetFieldPts(iso);
187  }
188 }
189 
190 
191 /*********************/
192 /* Function TwoPairs */
193 /*********************/
194 
198  int &pr)
199 /* The logic of the following SWITCH statement may only be
200  understood with a "peusdo truth-table" */
201 {
202  if (((cx[0]-cx[1])==0.0)&&
203  ((cy[0]-cy[1])==0.0)&&
204  ((cz[0]-cz[1])==0.0))
205  {
206  if (((cx[2]-cx[3])==0.0)&&
207  ((cy[2]-cy[3])==0.0)&&
208  ((cz[2]-cz[3])==0.0))
209  {
210  pr=4;
211  }
212  else
213  {
214  pr=3;
215  }
216  }
217  else
218  {
219  pr=1;
220  }
221 }
222 
223 
224 /*************************/
225 /* Function ThreeSimilar */
226 /*************************/
227 void ThreeSimilar (const int i, const int j, const int k,
228  int &pr, int &ps)
229 /* The logic of the following SWITCH statement may be
230  understood with a "peusdo truth-table" */
231 {
232  switch (i + j + k)
233  {
234  case (3):
235  pr = 3;
236  ps = 4;
237  break;
238  case (4):
239  pr = 2;
240  ps = 4;
241  break;
242  case (5):
243  if (j == 1)
244  {
245  pr = 2;
246  ps = 3;
247  }
248  else
249  {
250  pr = 1;
251  ps = 4;
252  }
253  break;
254  case (6):
255  if (i == 0)
256  {
257  pr = 1;
258  ps = 3;
259  }
260  else
261  {
262  pr = 0;
263  ps = 4;
264  }
265  break;
266  case (7):
267  if (i == 0)
268  {
269  pr = 1;
270  ps = 2;
271  }
272  else
273  {
274  pr = 0;
275  ps = 3;
276  }
277  break;
278  case (8):
279  pr = 0;
280  ps = 2;
281  break;
282  case (9):
283  pr = 0;
284  ps = 1;
285  break;
286  default:
287  ASSERTL0(false,"Error in 5-point triangulation in ThreeSimilar");
288  break;
289  }
290 }
291 
293  const int fieldid,
294  const NekDouble val)
295 {
296  vector<IsoSharedPtr> returnval;
297 
298  int coordim = m_f->m_exp[0]->GetCoordim(0);
299  int nfields = m_f->m_fieldPts->GetNFields() + coordim;
300 
301  ASSERTL0(coordim == 3,
302  "This methods is currently only set up for 3D fields");
303  ASSERTL1(coordim + fieldid < nfields,
304  "field id is larger than number contained in FieldPts");
306  m_f->m_fieldPts->GetPts(fields);
307 
308  Array<OneD, NekDouble> c = fields[coordim + fieldid];
309 
310  int i, j, k, ii, jj, kk, r, s, n, counter, boolean;
311  Array<OneD, Array<OneD, NekDouble> > intfields(nfields);
312  intfields[0] = Array<OneD, NekDouble>(5*nfields);
313  for(i = 1; i < nfields; ++i)
314  {
315  intfields[i] = intfields[i-1] + 5;
316  }
317  Array<OneD, NekDouble> cx = intfields[0];
318  Array<OneD, NekDouble> cy = intfields[1];
319  Array<OneD, NekDouble> cz = intfields[2];
320 
321  vector<Array<OneD, int> > ptsConn;
322  m_f->m_fieldPts->GetConnectivity(ptsConn);
323  for(int zone = 0; zone < ptsConn.size(); ++zone)
324  {
325  IsoSharedPtr iso;
326 
327  iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-3);
328 
329  int nelmt = ptsConn[zone].num_elements()
330  /(coordim+1);
331 
332  Array<OneD, int> conn = ptsConn[zone];
333 
334  for (n = 0, i = 0; i < nelmt; ++i)
335  {
336  // check to see if val is between vertex values
337  if(!(((c[conn[i*4]] >val)&&(c[conn[i*4+1]]>val)&&
338  (c[conn[i*4+2]]>val)&&(c[conn[i*4+3]]>val))||
339  ((c[conn[i*4 ]]<val)&&(c[conn[i*4+1]]<val)&&
340  (c[conn[i*4+2]]<val)&&(c[conn[i*4+3]]<val))))
341  {
342 
343  // loop over all edges and interpolate if
344  // contour is between vertex values
345  for (counter = 0, j=0; j<=2; j++)
346  {
347  for (k=j+1; k<=3; k++)
348  {
349  if (((c[conn[i*4+j]]>=val)&&
350  (val>=c[conn[i*4+k]]))||
351  ((c[conn[i*4+j]]<=val)&&
352  (val<=c[conn[i*4+k]])))
353  {
354  // linear interpolation of fields
355  // (and coords).
356  NekDouble cj = c[conn[i*4+j]];
357  NekDouble ck = c[conn[i*4+k]];
358  NekDouble factor = (val-cj)/(ck-cj);
359 
360  if(fabs(cj-ck) > 1e-12)
361  {
362  // interpolate coordinates and fields
363  for(int f = 0; f < nfields; ++f)
364  {
365  if(counter == 5)
366  {
367  ASSERTL0(false,"Counter is 5");
368  }
369  intfields[f][counter] =
370  fields[f][conn[4*i+j]] +
371  factor*(fields[f][conn[4*i+k]] -
372  fields[f][conn[4*i+j]]);
373  }
374  ++counter;
375  }
376  }
377  }
378  }
379 
380  switch(counter)
381  {
382  case 3:
383  n+=1;
384  iso->resize_fields(3*n);
385 
386  for(j = 0; j < 3; ++j)
387  {
388  iso->set_fields(3*(n-1)+j,intfields,j);
389  }
390  break;
391  case 4:
392  n+=2;
393  iso->resize_fields(3*n);
394 
395  for(j = 0; j < 3; ++j)
396  {
397  iso->set_fields(3*(n-2)+j,intfields,j);
398  iso->set_fields(3*(n-1)+j,intfields,j+1);
399  }
400  break;
401  case 5:
402  n+=1;
403  iso->resize_fields(3*n);
404 
405  boolean=0;
406  for (ii=0;ii<=2;ii++)
407  {
408  for (jj=ii+1;jj<=3;jj++)
409  {
410  for (kk=jj+1;kk<=4;kk++)
411  {
412  if((((cx[ii]-cx[jj])==0.0)&&
413  ((cy[ii]-cy[jj])==0.0)&&
414  ((cz[ii]-cz[jj])==0.0))&&
415  (((cx[ii]-cx[kk])==0.0)&&
416  ((cy[ii]-cy[kk])==0.0)&&
417  ((cz[ii]-cz[kk])==0.0)))
418  {
419  boolean+=1;
420  ThreeSimilar (ii,jj,kk,r,s);
421 
422  iso->set_fields(3*(n-1) ,intfields,ii);
423  iso->set_fields(3*(n-1)+1,intfields,r);
424  iso->set_fields(3*(n-1)+2,intfields,s);
425  }
426  else
427  {
428  boolean+=0;
429  }
430  }
431  }
432  }
433 
434  if (boolean==0)
435  {
436  TwoPairs (cx,cy,cz,r);
437 
438  iso->set_fields(3*(n-1) ,intfields,0);
439  iso->set_fields(3*(n-1)+1,intfields,2);
440  iso->set_fields(3*(n-1)+2,intfields,r);
441  }
442  break;
443  }
444  }
445  }
446  iso->set_ntris(n);
447 
448  // condense the information in this elemental extraction.
449  iso->condense();
450  returnval.push_back(iso);
451  }
452 
453  return returnval;
454 }
455 
456 // reset m_fieldPts with values from iso;
457 void ProcessIsoContour::ResetFieldPts(vector<IsoSharedPtr> &iso)
458 {
459  int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
460 
461  // set output to triangle block.
462  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
463 
464  Array<OneD, Array<OneD, NekDouble> > newfields(nfields);
465 
466  // count up number of points
467  int npts = 0;
468  for(int i =0; i < iso.size(); ++i)
469  {
470  npts += iso[i]->get_nvert();
471  }
472 
473  // set up coordinate in new field
474  newfields[0] = Array<OneD, NekDouble>(npts);
475  newfields[1] = Array<OneD, NekDouble>(npts);
476  newfields[2] = Array<OneD, NekDouble>(npts);
477 
478  int cnt = 0;
479  for(int i =0; i < iso.size(); ++i)
480  {
481  for(int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
482  {
483  newfields[0][cnt] = iso[i]->get_x(j);
484  newfields[1][cnt] = iso[i]->get_y(j);
485  newfields[2][cnt] = iso[i]->get_z(j);
486  }
487  }
488 
489 
490  // set up fields
491  for(int f = 0; f < nfields-3; ++f)
492  {
493  newfields[f+3] = Array<OneD, NekDouble>(npts);
494 
495  cnt = 0;
496  for(int i =0; i < iso.size(); ++i)
497  {
498  for(int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
499  {
500  newfields[f+3][cnt] = iso[i]->get_fields(f,j);
501  }
502  }
503  }
504 
505  m_f->m_fieldPts->SetPts(newfields);
506 
507  // set up connectivity data.
508  vector<Array<OneD, int> > ptsConn;
509  m_f->m_fieldPts->GetConnectivity(ptsConn);
510  cnt = 0;
511  ptsConn.clear();
512  for(int i =0; i < iso.size(); ++i)
513  {
514  int ntris = iso[i]->get_ntris();
515  Array<OneD, int> conn(ntris*3);
516 
517  for(int j = 0; j < 3*ntris; ++j)
518  {
519  conn[j] = cnt + iso[i]->get_vid(j);
520  }
521  ptsConn.push_back(conn);
522  cnt += iso[i]->get_nvert();
523  }
524  m_f->m_fieldPts->SetConnectivity(ptsConn);
525 }
526 
527 void Iso::condense(void)
528 {
529  register int i,j,cnt;
530  IsoVertex v;
531  vector<IsoVertex> vert;
533 
534  if(!m_ntris) return;
535 
536  if(m_condensed) return;
537  m_condensed = true;
538 
539  vert.reserve(m_ntris/6);
540 
542 
543  // fill first 3 points and initialise fields
544  v.m_fields.resize(m_fields.size());
545  for(cnt =0, i = 0; i < 3; ++i)
546  {
547  v.m_x = m_x[i];
548  v.m_y = m_y[i];
549  v.m_z = m_z[i];
550  for(int f = 0; f < m_fields.size(); ++f)
551  {
552  v.m_fields[f] = m_fields[f][i];
553  }
554  v.m_id = cnt;
555  vert.push_back(v);
556  m_vid[i] = v.m_id;
557  ++cnt;
558  }
559 
560  for(i = 1; i < m_ntris; ++i)
561  {
562  for(j = 0; j < 3; ++j)
563  {
564  v.m_x = m_x[3*i+j];
565  v.m_y = m_y[3*i+j];
566  v.m_z = m_z[3*i+j];
567 
568  pt = find(vert.begin(),vert.end(),v);
569  if(pt != vert.end())
570  {
571  m_vid[3*i+j] = pt[0].m_id;
572  }
573  else
574  {
575  v.m_id = cnt;
576 
577  for(int f = 0; f < m_fields.size(); ++f)
578  {
579  v.m_fields[f] = m_fields[f][3*i+j];
580  }
581 
582  vert.push_back(v);
583 
584  m_vid[3*i+j] = v.m_id;
585  ++cnt;
586  }
587  }
588  }
589 
590  // remove elements with multiple vertices
591  for(i = 0,cnt=0; i < m_ntris;)
592  {
593  if((m_vid[3*i] ==m_vid[3*i+1])||
594  (m_vid[3*i] ==m_vid[3*i+2])||
595  (m_vid[3*i+1]==m_vid[3*i+2]))
596  {
597  cnt++;
598  for(j = 3*i; j < 3*(m_ntris-1); ++j)
599  {
600  m_vid[j] = m_vid[j+3];
601  }
602  m_ntris--;
603  }
604  else
605  {
606  ++i;
607  }
608  }
609 
610  m_nvert = vert.size();
611 
612  m_x.resize(m_nvert);
613  m_y.resize(m_nvert);
614  m_z.resize(m_nvert);
615 
616  for(int f = 0; f < m_fields.size(); ++f)
617  {
618  m_fields[f].resize(m_nvert);
619  }
620 
621  for(i = 0; i < m_nvert; ++i)
622  {
623  m_x[i] = vert[i].m_x;
624  m_y[i] = vert[i].m_y;
625  m_z[i] = vert[i].m_z;
626  for(int f = 0; f < m_fields.size(); ++f)
627  {
628  m_fields[f][i] = vert[i].m_fields[f];
629  }
630  }
631 }
632 
633 
635 
636 // define == if point is within 1e-4
637 bool operator == (const IsoVertex& x, const IsoVertex& y)
638 {
639  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) +
640  (x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? true:false;
641 }
642 
643 // define != if point is outside 1e-4
644 bool operator != (const IsoVertex& x, const IsoVertex& y)
645 {
646  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) +
647  (x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? 0:1;
648 }
649 
650 
652  NekDouble x2, NekDouble y2, NekDouble z2)
653 {
654  if((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2) < SQ_PNT_TOL)
655  {
656  return true;
657  }
658 
659  return false;
660 }
661 
662 void Iso::globalcondense(vector<IsoSharedPtr> &iso)
663 {
664  int i,j,n;
665  int nvert,nelmt;
666  int niso=iso.size();
667  int id1,id2;
669 
670  if(m_condensed) return;
671  m_condensed = true;
672 
673  vidmap = Array<OneD, Array<OneD, int> > (niso);
674 
675  m_ntris = 0;
676  for(i = 0; i < niso; ++i)
677  {
678  if(iso[i]->m_ntris)
679  {
680  m_ntris += iso[i]->m_ntris;
681  }
682  }
683 
685 
686  m_nvert = 0;
687  for(i = 0; i < niso; ++i)
688  {
689  if(iso[i]->m_ntris)
690  {
691  m_nvert += iso[i]->m_nvert;
692  }
693  }
694 
695  vector< vector<int> > isocon;
696  isocon.resize(niso);
697 
698  // identify which iso are connected by at least one point;
699  // find min x,y,z and max x,y,z and see if overlap to select
700  // which zones should be connected
701  {
702  vector<Array<OneD, NekDouble> > sph(niso);
703  Array<OneD, NekDouble> rng(6);
704  for(i = 0; i < niso; ++i)
705  {
706  sph[i] = Array<OneD, NekDouble>(4);
707 
708  // find max and min of isocontour
709  rng[0] = rng[3] = iso[i]->m_x[0];
710  rng[1] = rng[4] = iso[i]->m_x[1];
711  rng[2] = rng[5] = iso[i]->m_x[2];
712 
713  for(id1 = 1; id1 < iso[i]->m_nvert;++id1)
714  {
715  rng[0] = min(rng[0],iso[i]->m_x[i]);
716  rng[1] = min(rng[1],iso[i]->m_y[i]);
717  rng[2] = min(rng[2],iso[i]->m_z[i]);
718 
719  rng[3] = max(rng[3],iso[i]->m_x[i]);
720  rng[4] = max(rng[4],iso[i]->m_y[i]);
721  rng[5] = max(rng[5],iso[i]->m_z[i]);
722  }
723 
724  // centroid
725  sph[i][0] = (rng[3]+rng[0])/2.0;
726  sph[i][1] = (rng[4]+rng[1])/2.0;
727  sph[i][2] = (rng[5]+rng[2])/2.0;
728 
729  // radius;
730  sph[i][3] = sqrt((rng[3]-rng[0])*(rng[3]-rng[0]) +
731  (rng[4]-rng[1])*(rng[4]-rng[1]) +
732  (rng[5]-rng[2])*(rng[5]-rng[2]));
733  }
734 
735  for(i = 0; i < niso; ++i)
736  {
737  for(j = i+1; j < niso; ++j)
738  {
739  NekDouble diff=sqrt((sph[i][0]-sph[j][0])*(sph[i][0]-sph[j][0])+
740  (sph[i][1]-sph[j][1])*(sph[i][1]-sph[j][1])+
741  (sph[i][2]-sph[j][2])*(sph[i][2]-sph[j][2]));
742 
743  // if centroid is closer than added radii
744  if(diff < sph[i][3] + sph[j][3])
745  {
746  isocon[i].push_back(j);
747  }
748  }
749  }
750  }
751 
752 
753  for(i = 0; i < niso; ++i)
754  {
755  vidmap[i] = Array<OneD, int>(iso[i]->m_nvert,-1);
756  }
757  nvert = 0;
758  // identify which vertices are connected to tolerance
759  cout << "GlobalCondense: Matching Vertices [" << endl << flush;
760  int cnt = 0;
761  // count up amount of checking to be done
762  NekDouble totpts = 0;
763  for(i = 0; i < niso; ++i)
764  {
765  for(n = 0; n < isocon[i].size(); ++n)
766  {
767  int con = isocon[i][n];
768  totpts += iso[i]->m_nvert*iso[con]->m_nvert;
769  }
770  }
771  int totchk = totpts/50;
772  int cnt_out = 0;
773  for(i = 0; i < niso; ++i)
774  {
775  for(n = 0; n < isocon[i].size(); ++n)
776  {
777  int con = isocon[i][n];
778  for(id1 = 0; id1 < iso[i]->m_nvert; ++id1)
779  {
780  for(id2 = 0; id2 < iso[con]->m_nvert; ++id2, ++cnt)
781  {
782  if(cnt%totchk == 0)
783  {
784  cout <<cnt_out << "%" << '\r' << flush;
785  cnt_out +=2;
786  }
787 
788  if((vidmap[con][id2] == -1)||(vidmap[i][id1] == -1))
789  {
790  if(same(iso[i]->m_x[id1], iso[i]->m_y[id1],
791  iso[i]->m_z[id1], iso[con]->m_x[id2],
792  iso[con]->m_y[id2],iso[con]->m_z[id2]))
793  {
794  if((vidmap[i][id1] == -1) &&
795  (vidmap[con][id2] != -1))
796  {
797  vidmap[i][id1] = vidmap[con][id2];
798  }
799  else if((vidmap[con][id2] == -1) &&
800  (vidmap[i][id1] != -1))
801  {
802  vidmap[con][id2] = vidmap[i][id1];
803  }
804  else if((vidmap[con][id2] == -1) &&
805  (vidmap[i][id1] == -1))
806  {
807  vidmap[i][id1] = vidmap[con][id2] = nvert++;
808  }
809  }
810  }
811  }
812  }
813  }
814 
815  for(id1 = 0; id1 < iso[i]->m_nvert;++id1)
816  {
817  if(vidmap[i][id1] == -1)
818  {
819  vidmap[i][id1] = nvert++;
820  }
821  }
822  }
823  cout <<endl << "]"<<endl;
824  m_nvert = nvert;
825 
826  nelmt = 0;
827  // reset m_vid;
828  for(n = 0; n < niso; ++n)
829  {
830  for(i = 0; i < iso[n]->m_ntris; ++i,nelmt++)
831  {
832  for(j=0; j < 3;++j)
833  {
834  m_vid[3*nelmt+j] = vidmap[n][iso[n]->m_vid[3*i+j]];
835  }
836  }
837  }
838 
839  m_ntris = nelmt;
840 
841  m_x.resize(m_nvert);
842  m_y.resize(m_nvert);
843  m_z.resize(m_nvert);
844 
845  m_fields.resize(iso[0]->m_fields.size());
846  for(i = 0; i < iso[0]->m_fields.size(); ++i)
847  {
848  m_fields[i].resize(m_nvert);
849  }
850 
851 
852  // reset coordinate and fields.
853  for(n = 0; n < niso; ++n)
854  {
855  for(i = 0; i < iso[n]->m_nvert; ++i)
856  {
857  m_x[vidmap[n][i]] = iso[n]->m_x[i];
858  m_y[vidmap[n][i]] = iso[n]->m_y[i];
859  m_z[vidmap[n][i]] = iso[n]->m_z[i];
860 
861  for(j = 0; j < m_fields.size(); ++j)
862  {
863  m_fields[j][vidmap[n][i]] = iso[n]->m_fields[j][i];
864  }
865  }
866  }
867 }
868 
869 void Iso::smooth(int n_iter, NekDouble lambda, NekDouble mu)
870 {
871  int iter,i,j;
872  NekDouble del_v[3];
873  NekDouble w;
874  Array<OneD, NekDouble> xtemp, ytemp, ztemp;
875  vector< vector<int> > adj,vertcon;
878 
879  // determine elements around each vertex
880  vertcon.resize(m_nvert);
881  for(i = 0; i < m_ntris; ++i)
882  {
883  for(j = 0; j < 3; ++j)
884  {
885  vertcon[m_vid[3*i+j]].push_back(i);
886  }
887  }
888 
889  // determine vertices around each vertex
890  adj.resize(m_nvert);
891 
892  for(i =0; i < m_nvert; ++i)
893  {
894  for(ipt = vertcon[i].begin(); ipt != vertcon[i].end(); ++ipt)
895  {
896  for(j = 0; j < 3; ++j)
897  {
898  // make sure not adding own vertex
899  if(m_vid[3*(*ipt)+j] != i)
900  {
901  // check to see if vertex has already been added
902  for(iad = adj[i].begin(); iad != adj[i].end();++iad)
903  {
904  if(*iad == (m_vid[3*(*ipt)+j])) break;
905  }
906 
907  if(iad == adj[i].end())
908  {
909  adj[i].push_back(m_vid[3*(*ipt)+j]);
910  }
911  }
912  }
913  }
914  }
915 
919 
920  // smooth each point
921  for (iter=0;iter<n_iter;iter++)
922  {
923  // compute first weighted average
924  for(i=0;i< m_nvert;++i)
925  {
926  w = 1.0/(NekDouble)(adj[i].size());
927 
928  del_v[0] = del_v[1] = del_v[2] = 0.0;
929 
930  for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
931  {
932  del_v[0] = del_v[0] + (m_x[*iad]-m_x[i])*w;
933  del_v[1] = del_v[1] + (m_y[*iad]-m_y[i])*w;
934  del_v[2] = del_v[2] + (m_z[*iad]-m_z[i])*w;
935  }
936 
937  xtemp[i] = m_x[i] + del_v[0] * lambda;
938  ytemp[i] = m_y[i] + del_v[1] * lambda;
939  ztemp[i] = m_z[i] + del_v[2] * lambda;
940  }
941 
942  // compute second weighted average
943  for(i=0;i< m_nvert;++i)
944  {
945 
946  w = 1.0/(NekDouble)(adj[i].size());
947  del_v[0] = del_v[1] = del_v[2] = 0;
948 
949  for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
950  {
951  del_v[0] = del_v[0] + (m_x[*iad]-m_x[i])*w;
952  del_v[1] = del_v[1] + (m_y[*iad]-m_y[i])*w;
953  del_v[2] = del_v[2] + (m_z[*iad]-m_z[i])*w;
954  }
955 
956  m_x[i] = xtemp[i] + del_v[0] * mu;
957  m_y[i] = ytemp[i] + del_v[1] * mu;
958  m_z[i] = ztemp[i] + del_v[2] * mu;
959  }
960  }
961 }
962 
963 
964  void Iso::separate_regions(vector<IsoSharedPtr> &sep_iso, int minsize)
965  {
966  int i,j,k,id;
969  list<int> tocheck;
971 
972  Array<OneD, bool> viddone(m_nvert,false);
973 
974  // make list of connecting tris around each vertex
975  for(i = 0; i < m_ntris; ++i)
976  {
977  for(j = 0; j < 3; ++j)
978  {
979  vertcon[m_vid[3*i+j]].push_back(i);
980  }
981  }
982 
983  Array<OneD, int> vidregion(m_nvert,-1);
984 
985  int nregions = -1;
986 
987  cout << "Identifying separate regions [." << flush ;
988  // check all points are assigned;
989  for(k = 0; k < m_nvert; ++k)
990  {
991  if(vidregion[k] == -1)
992  {
993  vidregion[k] = ++nregions;
994 
995  // find all elmts around this.. vertex that need to be checked
996  for(ipt = vertcon[k].begin(); ipt != vertcon[k].end(); ++ipt)
997  {
998  for(i = 0; i < 3; ++i)
999  {
1000  if(vidregion[id = m_vid[3*(ipt[0])+i]] == -1)
1001  {
1002  tocheck.push_back(id);
1003  vidregion[id] = nregions;
1004  }
1005  }
1006  }
1007  viddone[k] = 1;
1008 
1009  // check all other neighbouring vertices
1010  while(tocheck.size())
1011  {
1012  cid = tocheck.begin();
1013  while(cid != tocheck.end())
1014  {
1015  if(!viddone[*cid])
1016  {
1017  for(ipt = vertcon[*cid].begin(); ipt != vertcon[*cid].end(); ++ipt)
1018  {
1019  for(i = 0; i < 3; ++i)
1020  {
1021  if(vidregion[id = m_vid[3*(ipt[0])+i]] == -1)
1022  {
1023  tocheck.push_back(id);
1024  vidregion[id] = nregions;
1025  }
1026  }
1027  }
1028  viddone[*cid] = 1;
1029  ++cid;
1030  tocheck.pop_front();
1031  }
1032  }
1033  }
1034  }
1035  }
1036  nregions++;
1037 
1038  cout << "]" << endl << flush ;
1039 
1040  Array<OneD, int> nvert(nregions,0);
1041  Array<OneD, int> nelmt(nregions,0);
1042 
1043  // count nverts
1044  for(i = 0; i < m_nvert; ++i)
1045  {
1046  nvert[vidregion[i]] +=1;
1047  }
1048 
1049  // count nelmts
1050  for(i = 0; i < m_ntris; ++i)
1051  {
1052  nelmt[vidregion[m_vid[3*i]]] +=1;
1053  }
1054 
1055  Array<OneD, int> vidmap(m_nvert);
1056  // generate new list of isocontour
1057  for(int n = 0; n < nregions; ++n)
1058  {
1059  if(nelmt[n] > minsize)
1060  {
1061  int nfields = m_fields.size();
1063 
1064  iso->m_ntris = nelmt[n];
1065  iso->m_vid = Array<OneD, int>(3*nelmt[n]);
1066 
1067  iso->m_nvert = nvert[n];
1068  iso->m_x.resize(nvert[n]);
1069  iso->m_y.resize(nvert[n]);
1070  iso->m_z.resize(nvert[n]);
1071 
1072  iso->m_fields.resize(nfields);
1073  for(i = 0; i < nfields; ++i)
1074  {
1075  iso->m_fields[i].resize(nvert[n]);
1076  }
1077 
1078 
1079  int cnt = 0;
1080  // generate vid map;
1081  Vmath::Zero(m_nvert,vidmap,1);
1082  for(i = 0; i < m_nvert; ++i)
1083  {
1084  if(vidregion[i] == n)
1085  {
1086  vidmap[i] = cnt++;
1087  }
1088  }
1089 
1090  cnt = 0;
1091  for(i = 0; i < m_ntris; ++i)
1092  {
1093  if(vidregion[m_vid[3*i]] == n)
1094  {
1095  for(j = 0; j < 3; ++j)
1096  {
1097  iso->m_vid[3*cnt+j] = vidmap[m_vid[3*i+j]];
1098 
1099  iso->m_x[vidmap[m_vid[3*i+j]]] = m_x[m_vid[3*i+j]];
1100  iso->m_y[vidmap[m_vid[3*i+j]]] = m_y[m_vid[3*i+j]];
1101  iso->m_z[vidmap[m_vid[3*i+j]]] = m_z[m_vid[3*i+j]];
1102 
1103  for(k = 0; k < nfields; ++k)
1104  {
1105  iso->m_fields[k][vidmap[m_vid[3*i+j]]] = m_fields[k][m_vid[3*i+j]];
1106  }
1107  }
1108  cnt++;
1109  }
1110  }
1111 
1112  WARNINGL0(cnt == nelmt[n],"Number of elements do not match");
1113  sep_iso.push_back(iso);
1114  }
1115  }
1116  }
1117 }
1118 }
1119 
1120 
#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.
void globalcondense(vector< boost::shared_ptr< Iso > > &iso)
FieldSharedPtr m_f
Field object.
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 separate_regions(vector< boost::shared_ptr< Iso > > &iso, int minsize)
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:677
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
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:312
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.
Definition: Vmath.cpp:359
#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.
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215