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: Set up fields as interpolation to equispaced output
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 #include <string>
36 #include <iostream>
37 using namespace std;
38 
39 #include "ProcessIsoContour.h"
40 
41 #include <LibUtilities/BasicUtils/SharedArray.hpp>
42 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
43 #include <boost/math/special_functions/fpclassify.hpp>
44 
45 namespace Nektar
46 {
47 namespace Utilities
48 {
49 
50 ModuleKey ProcessIsoContour::className =
51  GetModuleFactory().RegisterCreatorFunction(
52  ModuleKey(eProcessModule, "isocontour"),
53  ProcessIsoContour::create,
54  "Extract an isocontour of fieldid variable and at value fieldvalue, Optionally fieldstr can be specified for a string defiition or smooth for smoothing");
55 
56 ProcessIsoContour::ProcessIsoContour(FieldSharedPtr f) :
57  ProcessEquiSpacedOutput(f)
58 {
59 
60  m_config["fieldstr"] = ConfigOption(false, "NotSet",
61  "string of isocontour to be extracted");
62  m_config["fieldname"] = ConfigOption(false, "isocon",
63  "name for isocontour if fieldstr "
64  "specified, default is isocon");
65 
66  m_config["fieldid"] = ConfigOption(false, "NotSet",
67  "field id to extract");
68  m_config["fieldvalue"] = ConfigOption(false, "NotSet",
69  "field value to extract");
70 
71  m_config["globalcondense"] = ConfigOption(true, "NotSet",
72  "Globally condense contour to unique "
73  "values");
74 
75  m_config["smooth"] = ConfigOption(true, "NotSet",
76  "Smooth isocontour (implies global "
77  "condense)");
78  m_config["smoothiter"] = ConfigOption(false, "100",
79  "Number of smoothing cycle, default = "
80  "100");
81  m_config["smoothposdiffusion"] = ConfigOption(false,"0.5",
82  "Postive diffusion coefficient "
83  "(0 < lambda < 1), default = 0.5");
84  m_config["smoothnegdiffusion"] = ConfigOption(false,"0.495",
85  "Negative diffusion coefficient "
86  "(0 < mu < 1), default = 0.495");
87 
88 }
89 
90 ProcessIsoContour::~ProcessIsoContour(void)
91 {
92 }
93 
94 void ProcessIsoContour::Process(po::variables_map &vm)
95 {
96  vector<IsoSharedPtr> iso;
97 
98  // extract all fields to equi-spaced
99  SetupEquiSpacedField();
100 
101  int fieldid;
102  NekDouble value;
103 
104  if(m_config["fieldstr"].m_beenSet) //generate field of interest
105  {
106 
107  m_f->m_fieldPts->m_nFields++;
108 
109  int nfields = m_f->m_fieldPts->m_nFields;
110  int coordim = m_f->m_fieldPts->m_ptsDim;
111 
112  fieldid = nfields-1;
113 
114  Array<OneD, Array<OneD, NekDouble> > newpts(nfields+coordim);
115 
116  // redirect existing pts
117  for(int i = 0; i < nfields+coordim-1; ++i)
118  {
119  newpts[i] = m_f->m_fieldPts->m_pts[i];
120  }
121  int npts = m_f->m_fieldPts->m_pts[0].num_elements();
122  newpts[nfields+coordim-1] = Array<OneD, NekDouble>(npts);
123 
124  m_f->m_fieldPts->m_pts = newpts;
125 
126  // evaluate new function
127  LibUtilities::AnalyticExpressionEvaluator strEval;
128  string varstr = "x y z";
129  vector<Array<OneD, const NekDouble> > interpfields;
130 
131  for(int i = 0; i < coordim; ++i)
132  {
133  interpfields.push_back(m_f->m_fieldPts->m_pts[i]);
134  }
135  for(int i = 0; i < nfields-1; ++i)
136  {
137  varstr += " " + m_f->m_fieldPts->m_fields[i];
138  interpfields.push_back(m_f->m_fieldPts->m_pts[i+3]);
139  }
140 
141  int ExprId = -1;
142  std::string str = m_config["fieldstr"].as<string>();
143  ExprId = strEval.DefineFunction(varstr.c_str(), str);
144 
145  strEval.Evaluate(ExprId, interpfields,
146  m_f->m_fieldPts->m_pts[nfields+coordim-1]);
147 
148  // set up field name if provided otherwise called "isocon" from default
149  m_f->m_fieldPts->m_fields.push_back(
150  m_config["fieldname"].as<string>());
151  }
152  else
153  {
154  fieldid = m_config["fieldid"].as<int>();
155  }
156 
157  value = m_config["fieldvalue"].as<NekDouble>();
158 
159  iso = ExtractContour(fieldid,value);
160  bool smoothing = m_config["smooth"].m_beenSet;
161  bool globalcondense = m_config["globalcondense"].m_beenSet;
162  if(smoothing||globalcondense)
163  {
164  vector<IsoSharedPtr> glob_iso;
165  int nfields = m_f->m_fieldPts->m_pts.num_elements();
166  IsoSharedPtr g_iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-3);
167 
168  g_iso->globalcondense(iso);
169 
170  if(smoothing)
171  {
172  int niter = m_config["smoothiter"].as<int>();
173  NekDouble lambda = m_config["smoothposdiffusion"].as<NekDouble>();
174  NekDouble mu = m_config["smoothnegdiffusion"].as<NekDouble>();
175  g_iso->smooth(niter,lambda,-mu);
176  }
177 
178  glob_iso.push_back(g_iso);
179 
180  ResetFieldPts(glob_iso);
181  }
182  else
183  {
184  ResetFieldPts(iso);
185  }
186 }
187 
188 
189 /*********************/
190 /* Function TwoPairs */
191 /*********************/
192 
193 void TwoPairs (Array<OneD, NekDouble> &cx,
194  Array<OneD, NekDouble> &cy,
195  Array<OneD, NekDouble> &cz,
196  int &pr)
197 /* The logic of the following SWITCH statement may only be
198  understood with a "peusdo truth-table" */
199 {
200  if (((cx[0]-cx[1])==0.0)&&
201  ((cy[0]-cy[1])==0.0)&&
202  ((cz[0]-cz[1])==0.0))
203  {
204  if (((cx[2]-cx[3])==0.0)&&
205  ((cy[2]-cy[3])==0.0)&&
206  ((cz[2]-cz[3])==0.0))
207  {
208  pr=4;
209  } else
210  {
211  pr=3;
212  }
213  }
214  else
215  {
216  pr=1;
217  }
218 }
219 
220 
221 /*************************/
222 /* Function ThreeSimilar */
223 /*************************/
224 void ThreeSimilar (const int i, const int j, const int k,
225  int &pr, int &ps)
226 /* The logic of the following SWITCH statement may be
227  understood with a "peusdo truth-table" */
228 {
229  switch (i + j + k)
230  {
231  case (3):
232  pr = 3;
233  ps = 4;
234  break;
235  case (4):
236  pr = 2;
237  ps = 4;
238  break;
239  case (5):
240  if (j == 1)
241  {
242  pr = 2;
243  ps = 3;
244  }
245  else
246  {
247  pr = 1;
248  ps = 4;
249  }
250  break;
251  case (6):
252  if (i == 0)
253  {
254  pr = 1;
255  ps = 3;
256  }
257  else
258  {
259  pr = 0;
260  ps = 4;
261  }
262  break;
263  case (7):
264  if (i == 0)
265  {
266  pr = 1;
267  ps = 2;
268  }
269  else
270  {
271  pr = 0;
272  ps = 3;
273  }
274  break;
275  case (8):
276  pr = 0;
277  ps = 2;
278  break;
279  case (9):
280  pr = 0;
281  ps = 1;
282  break;
283  default:
284  printf("Error in 5-point triangulation in ThreeSimilar");
285  break;
286  }
287 }
288 
289 vector<IsoSharedPtr> ProcessIsoContour::ExtractContour(
290  const int fieldid,
291  const NekDouble val)
292 {
293  vector<IsoSharedPtr> returnval;
294 
295  Array<OneD, NekDouble> x = m_f->m_fieldPts->m_pts[0];
296  Array<OneD, NekDouble> y = m_f->m_fieldPts->m_pts[1];
297  Array<OneD, NekDouble> z = m_f->m_fieldPts->m_pts[2];
298 
299  int coordim = m_f->m_exp[0]->GetCoordim(0);
300  int nfields = m_f->m_fieldPts->m_pts.num_elements();
301 
302  ASSERTL0(coordim == 3,
303  "This methods is currently only set up for 3D fields");
304  ASSERTL1(coordim + fieldid < nfields,
305  "field id is larger than number contained in FieldPts");
306  Array<OneD, Array<OneD, NekDouble> > fields = m_f->m_fieldPts->m_pts;
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  for(int zone = 0; zone < m_f->m_fieldPts->m_ptsConn.size(); ++zone)
322  {
323  IsoSharedPtr iso;
324 
325  iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-3);
326 
327  int nelmt = m_f->m_fieldPts->m_ptsConn[zone].num_elements()
328  /(coordim+1);
329 
330  Array<OneD, int> conn = m_f->m_fieldPts->m_ptsConn[zone];
331 
332  for (n = 0, i = 0; i < nelmt; ++i)
333  {
334  // check to see if val is between vertex values
335  if(!(((c[conn[i*4]] >val)&&(c[conn[i*4+1]]>val)&&
336  (c[conn[i*4+2]]>val)&&(c[conn[i*4+3]]>val))||
337  ((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  {
340 
341  // loop over all edges and interpolate if
342  // contour is between vertex values
343  for (counter = 0, j=0; j<=2; j++)
344  {
345  for (k=j+1; k<=3; k++)
346  {
347  if (((c[conn[i*4+j]]>=val)&&
348  (val>=c[conn[i*4+k]]))||
349  ((c[conn[i*4+j]]<=val)&&
350  (val<=c[conn[i*4+k]])))
351  {
352  // linear interpolation of fields
353  // (and coords).
354  NekDouble cj = c[conn[i*4+j]];
355  NekDouble ck = c[conn[i*4+k]];
356  NekDouble factor = (val-cj)/(ck-cj);
357 
358  if(fabs(cj-ck) > 1e-12)
359  {
360  // interpolate coordinates and fields
361  for(int f = 0; f < nfields; ++f)
362  {
363  if(counter == 5)
364  {
365  ASSERTL0(false,"Counter is 5");
366  }
367  intfields[f][counter] =
368  fields[f][conn[4*i+j]] +
369  factor*(fields[f][conn[4*i+k]] -
370  fields[f][conn[4*i+j]]);
371  }
372  ++counter;
373  }
374  }
375  }
376  }
377 
378  switch(counter)
379  {
380  case 3:
381  n+=1;
382  iso->resize_fields(3*n);
383 
384  for(j = 0; j < 3; ++j)
385  {
386  iso->set_fields(3*(n-1)+j,intfields,j);
387 
388  }
389  break;
390  case 4:
391  n+=2;
392  iso->resize_fields(3*n);
393 
394  for(j = 0; j < 3; ++j)
395  {
396  iso->set_fields(3*(n-2)+j,intfields,j);
397  iso->set_fields(3*(n-1)+j,intfields,j+1);
398 
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->m_pts.num_elements();
460 
461  // set output to triangle block.
462  m_f->m_fieldPts->m_ptype = 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->m_pts = newfields;
506 
507  // set up connectivity data.
508  cnt = 0;
509  m_f->m_fieldPts->m_ptsConn.clear();
510  for(int i =0; i < iso.size(); ++i)
511  {
512  int ntris = iso[i]->get_ntris();
513  Array<OneD, int> conn(ntris*3);
514 
515  for(int j = 0; j < 3*ntris; ++j)
516  {
517  conn[j] = cnt + iso[i]->get_vid(j);
518  }
519  m_f->m_fieldPts->m_ptsConn.push_back(conn);
520  cnt += iso[i]->get_nvert();
521  }
522 }
523 
524 void Iso::condense(void)
525 {
526  register int i,j,cnt;
527  IsoVertex v;
528  vector<IsoVertex> vert;
529  vector<IsoVertex>::iterator pt;
530 
531  if(!m_ntris) return;
532 
533  if(m_condensed) return;
534  m_condensed = true;
535 
536  vert.reserve(m_ntris/6);
537 
538  m_vid = Array<OneD, int>(3*m_ntris);
539 
540  // fill first 3 points and initialise fields
541  v.m_fields.resize(m_fields.size());
542  for(cnt =0, i = 0; i < 3; ++i)
543  {
544  v.m_x = m_x[i];
545  v.m_y = m_y[i];
546  v.m_z = m_z[i];
547  for(int f = 0; f < m_fields.size(); ++f)
548  {
549  v.m_fields[f] = m_fields[f][i];
550  }
551  v.m_id = cnt;
552  vert.push_back(v);
553  m_vid[i] = v.m_id;
554  ++cnt;
555  }
556 
557  for(i = 1; i < m_ntris; ++i)
558  {
559  for(j = 0; j < 3; ++j)
560  {
561  v.m_x = m_x[3*i+j];
562  v.m_y = m_y[3*i+j];
563  v.m_z = m_z[3*i+j];
564 
565  pt = find(vert.begin(),vert.end(),v);
566  if(pt != vert.end())
567  {
568  m_vid[3*i+j] = pt[0].m_id;
569  }
570  else
571  {
572  v.m_id = cnt;
573 
574  for(int f = 0; f < m_fields.size(); ++f)
575  {
576  v.m_fields[f] = m_fields[f][3*i+j];
577  }
578 
579  vert.push_back(v);
580 
581  m_vid[3*i+j] = v.m_id;
582  ++cnt;
583  }
584  }
585  }
586 
587  // remove elements with multiple vertices
588  for(i = 0,cnt=0; i < m_ntris;)
589  {
590  if((m_vid[3*i] ==m_vid[3*i+1])||
591  (m_vid[3*i] ==m_vid[3*i+2])||
592  (m_vid[3*i+1]==m_vid[3*i+2]))
593  {
594  cnt++;
595  for(j = 3*i; j < 3*(m_ntris-1); ++j)
596  {
597  m_vid[j] = m_vid[j+3];
598  }
599  m_ntris--;
600  }
601  else
602  {
603  ++i;
604  }
605  }
606 
607  m_nvert = vert.size();
608 
609  m_x.resize(m_nvert);
610  m_y.resize(m_nvert);
611  m_z.resize(m_nvert);
612 
613  for(int f = 0; f < m_fields.size(); ++f)
614  {
615  m_fields[f].resize(m_nvert);
616  }
617 
618  for(i = 0; i < m_nvert; ++i)
619  {
620  m_x[i] = vert[i].m_x;
621  m_y[i] = vert[i].m_y;
622  m_z[i] = vert[i].m_z;
623  for(int f = 0; f < m_fields.size(); ++f)
624  {
625  m_fields[f][i] = vert[i].m_fields[f];
626  }
627  }
628 }
629 
630 
631 NekDouble SQ_PNT_TOL=1e-16;
632 
633 // define == if point is within 1e-4
634 bool operator == (const IsoVertex& x, const IsoVertex& y)
635 {
636  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) +
637  (x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? true:false;
638 }
639 
640 // define != if point is outside 1e-4
641 bool operator != (const IsoVertex& x, const IsoVertex& y)
642 {
643  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) +
644  (x.m_z-y.m_z)*(x.m_z-y.m_z) < SQ_PNT_TOL)? 0:1;
645 }
646 
647 
648 bool same(NekDouble x1, NekDouble y1, NekDouble z1,
649  NekDouble x2, NekDouble y2, NekDouble z2)
650 {
651  if((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2) < SQ_PNT_TOL)
652  {
653  return true;
654  }
655 
656  return false;
657 }
658 
659 void Iso::globalcondense(vector<IsoSharedPtr> &iso)
660 {
661  int i,j,n;
662  int nvert,nelmt;
663  int niso=iso.size();
664  int id1,id2;
665  Array<OneD, Array<OneD, int> > vidmap;
666 
667  if(m_condensed) return;
668  m_condensed = true;
669 
670  vidmap = Array<OneD, Array<OneD, int> > (niso);
671 
672  m_ntris = 0;
673  for(i = 0; i < niso; ++i)
674  {
675  if(iso[i]->m_ntris)
676  {
677  m_ntris += iso[i]->m_ntris;
678  }
679  }
680 
681  m_vid = Array<OneD, int>(3*m_ntris);
682 
683  m_nvert = 0;
684  for(i = 0; i < niso; ++i)
685  {
686  if(iso[i]->m_ntris)
687  {
688  m_nvert += iso[i]->m_nvert;
689  }
690  }
691 
692  vector< vector<int> > isocon;
693  isocon.resize(niso);
694 
695  // identify which iso are connected by at least one point;
696  // find min x,y,z and max x,y,z and see if overlap to select
697  // which zones should be connected
698  {
699  vector<Array<OneD, NekDouble> > sph(niso);
700  Array<OneD, NekDouble> rng(6);
701  for(i = 0; i < niso; ++i)
702  {
703  sph[i] = Array<OneD, NekDouble>(4);
704 
705  // find max and min of isocontour
706  rng[0] = rng[3] = iso[i]->m_x[0];
707  rng[1] = rng[4] = iso[i]->m_x[1];
708  rng[2] = rng[5] = iso[i]->m_x[2];
709 
710  for(id1 = 1; id1 < iso[i]->m_nvert;++id1)
711  {
712  rng[0] = min(rng[0],iso[i]->m_x[i]);
713  rng[1] = min(rng[1],iso[i]->m_y[i]);
714  rng[2] = min(rng[2],iso[i]->m_z[i]);
715 
716  rng[3] = max(rng[3],iso[i]->m_x[i]);
717  rng[4] = max(rng[4],iso[i]->m_y[i]);
718  rng[5] = max(rng[5],iso[i]->m_z[i]);
719  }
720 
721  // centroid
722  sph[i][0] = (rng[3]+rng[0])/2.0;
723  sph[i][1] = (rng[4]+rng[1])/2.0;
724  sph[i][2] = (rng[5]+rng[2])/2.0;
725 
726  // radius;
727  sph[i][3] = sqrt((rng[3]-rng[0])*(rng[3]-rng[0]) +
728  (rng[4]-rng[1])*(rng[4]-rng[1]) +
729  (rng[5]-rng[2])*(rng[5]-rng[2]));
730  }
731 
732  for(i = 0; i < niso; ++i)
733  {
734  for(j = i+1; j < niso; ++j)
735  {
736  NekDouble diff=sqrt((sph[i][0]-sph[j][0])*(sph[i][0]-sph[j][0])+
737  (sph[i][1]-sph[j][1])*(sph[i][1]-sph[j][1])+
738  (sph[i][2]-sph[j][2])*(sph[i][2]-sph[j][2]));
739 
740  // if centroid is closer than added radii
741  if(diff < sph[i][3] + sph[j][3])
742  {
743  isocon[i].push_back(j);
744  }
745  }
746  }
747  }
748 
749 
750  for(i = 0; i < niso; ++i)
751  {
752  vidmap[i] = Array<OneD, int>(iso[i]->m_nvert,-1);
753  }
754  nvert = 0;
755  // identify which vertices are connected to tolerance
756  cout << "Matching Vertices [" <<flush;
757  int cnt = 0;
758  for(i = 0; i < niso; ++i)
759  {
760  for(n = 0; n < isocon[i].size(); ++n)
761  {
762  int con = isocon[i][n];
763  for(id1 = 0; id1 < iso[i]->m_nvert;++id1)
764  {
765  for(id2 = 0; id2 < iso[con]->m_nvert;++id2,++cnt)
766  {
767  if(cnt%1000000 == 0)
768  {
769  cout <<"." <<flush;
770  }
771 
772  //if((vidmap[con][id2] != -1)&&(vidmap[i][id1] != -1))
773  {
774  if(same(iso[i]->m_x[id1], iso[i]->m_y[id1],
775  iso[i]->m_z[id1], iso[con]->m_x[id2],
776  iso[con]->m_y[id2],iso[con]->m_z[id2]))
777  {
778  if((vidmap[i][id1] == -1) &&
779  (vidmap[con][id2] != -1))
780  {
781  vidmap[i][id1] = vidmap[con][id2];
782  }
783  else if((vidmap[con][id2] == -1) &&
784  (vidmap[i][id1] != -1))
785  {
786  vidmap[con][id2] = vidmap[i][id1];
787  }
788  else if((vidmap[con][id2] == -1) &&
789  (vidmap[i][id1] == -1))
790  {
791  vidmap[i][id1] = vidmap[con][id2] = nvert++;
792  }
793  }
794  }
795  }
796  }
797  }
798 
799  for(id1 = 0; id1 < iso[i]->m_nvert;++id1)
800  {
801  if(vidmap[i][id1] == -1)
802  {
803  vidmap[i][id1] = nvert++;
804  }
805  }
806  }
807  cout <<"]"<<endl;
808  m_nvert = nvert;
809 
810  nelmt = 0;
811  // reset m_vid;
812  for(n = 0; n < niso; ++n)
813  {
814  for(i = 0; i < iso[n]->m_ntris; ++i,nelmt++)
815  {
816  for(j=0; j < 3;++j)
817  {
818  m_vid[3*nelmt+j] = vidmap[n][iso[n]->m_vid[3*i+j]];
819  }
820  }
821  }
822 
823  m_ntris = nelmt;
824 
825  m_x.resize(m_nvert);
826  m_y.resize(m_nvert);
827  m_z.resize(m_nvert);
828 
829  m_fields.resize(iso[0]->m_fields.size());
830  for(i = 0; i < iso[0]->m_fields.size(); ++i)
831  {
832  m_fields[i].resize(m_nvert);
833  }
834 
835 
836  // reset coordinate and fields.
837  for(n = 0; n < niso; ++n)
838  {
839  for(i = 0; i < iso[n]->m_nvert; ++i)
840  {
841  m_x[vidmap[n][i]] = iso[n]->m_x[i];
842  m_y[vidmap[n][i]] = iso[n]->m_y[i];
843  m_z[vidmap[n][i]] = iso[n]->m_z[i];
844 
845  for(j = 0; j < m_fields.size(); ++j)
846  {
847  m_fields[j][vidmap[n][i]] = iso[n]->m_fields[j][i];
848  }
849  }
850  }
851 }
852 
853 void Iso::smooth(int n_iter, NekDouble lambda, NekDouble mu)
854 {
855  int iter,i,j;
856  NekDouble del_v[3];
857  NekDouble w;
858  Array<OneD, NekDouble> xtemp, ytemp, ztemp;
859  vector< vector<int> > adj,vertcon;
860  vector<int>::iterator iad;
861  vector<int>::iterator ipt;
862 
863  // determine elements around each vertex
864  vertcon.resize(m_nvert);
865  for(i = 0; i < m_ntris; ++i)
866  {
867  for(j = 0; j < 3; ++j)
868  {
869  vertcon[m_vid[3*i+j]].push_back(i);
870  }
871  }
872 
873  // determine vertices around each vertex
874  adj.resize(m_nvert);
875 
876  for(i =0; i < m_nvert; ++i)
877  {
878  for(ipt = vertcon[i].begin(); ipt != vertcon[i].end(); ++ipt)
879  {
880  for(j = 0; j < 3; ++j)
881  {
882  // make sure not adding own vertex
883  if(m_vid[3*(*ipt)+j] != i)
884  {
885  // check to see if vertex has already been added
886  for(iad = adj[i].begin(); iad != adj[i].end();++iad)
887  {
888  if(*iad == (m_vid[3*(*ipt)+j])) break;
889  }
890 
891  if(iad == adj[i].end())
892  {
893  adj[i].push_back(m_vid[3*(*ipt)+j]);
894  }
895  }
896  }
897  }
898  }
899 
900  xtemp = Array<OneD, NekDouble>(m_nvert);
901  ytemp = Array<OneD, NekDouble>(m_nvert);
902  ztemp = Array<OneD, NekDouble>(m_nvert);
903 
904  // smooth each point
905  for (iter=0;iter<n_iter;iter++)
906  {
907  // compute first weighted average
908  for(i=0;i< m_nvert;++i)
909  {
910  w = 1.0/(NekDouble)(adj[i].size());
911 
912  del_v[0] = del_v[1] = del_v[2] = 0.0;
913 
914  for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
915  {
916  del_v[0] = del_v[0] + (m_x[*iad]-m_x[i])*w;
917  del_v[1] = del_v[1] + (m_y[*iad]-m_y[i])*w;
918  del_v[2] = del_v[2] + (m_z[*iad]-m_z[i])*w;
919  }
920 
921  xtemp[i] = m_x[i] + del_v[0] * lambda;
922  ytemp[i] = m_y[i] + del_v[1] * lambda;
923  ztemp[i] = m_z[i] + del_v[2] * lambda;
924  }
925 
926  // compute second weighted average
927  for(i=0;i< m_nvert;++i)
928  {
929 
930  w = 1.0/(NekDouble)(adj[i].size());
931  del_v[0] = del_v[1] = del_v[2] = 0;
932 
933  for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
934  {
935  del_v[0] = del_v[0] + (m_x[*iad]-m_x[i])*w;
936  del_v[1] = del_v[1] + (m_y[*iad]-m_y[i])*w;
937  del_v[2] = del_v[2] + (m_z[*iad]-m_z[i])*w;
938  }
939 
940  m_x[i] = xtemp[i] + del_v[0] * mu;
941  m_y[i] = ytemp[i] + del_v[1] * mu;
942  m_z[i] = ztemp[i] + del_v[2] * mu;
943  }
944  }
945 }
946 
947 }
948 }
949 
950 
Generated on Wed Mar 4 2015 16:33:17 for Nektar++ by   doxygen 1.8.1.2