This processing module extracts an isocontour. More...

#include <ProcessIsoContour.h>

Inheritance diagram for Nektar::Utilities::ProcessIsoContour:

Collaboration diagram for Nektar::Utilities::ProcessIsoContour:

Public Member Functions
	ProcessIsoContour (FieldSharedPtr f)

virtual	~ProcessIsoContour ()

virtual void	Process (po::variables_map &vm)
	Write mesh to output file. More...

Public Member Functions inherited from Nektar::Utilities::ProcessEquiSpacedOutput
	ProcessEquiSpacedOutput (FieldSharedPtr f)

virtual	~ProcessEquiSpacedOutput ()

Public Member Functions inherited from Nektar::Utilities::ProcessModule
	ProcessModule ()

	ProcessModule (FieldSharedPtr p_f)

	ProcessModule (MeshSharedPtr p_m)

Public Member Functions inherited from Nektar::Utilities::Module
	Module (FieldSharedPtr p_f)

void	RegisterConfig (string key, string value)
	Register a configuration option with a module. More...

void	PrintConfig ()
	Print out all configuration options for a module. More...

void	SetDefaults ()
	Sets default configuration options for those which have not been set. More...

bool	GetRequireEquiSpaced (void)

void	SetRequireEquiSpaced (bool pVal)

void	EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)

	Module (MeshSharedPtr p_m)

virtual void	Process ()=0

void	RegisterConfig (std::string key, std::string value)

void	PrintConfig ()

void	SetDefaults ()

MeshSharedPtr	GetMesh ()

virtual void	ProcessVertices ()
	Extract element vertices. More...

virtual void	ProcessEdges (bool ReprocessEdges=true)
	Extract element edges. More...

virtual void	ProcessFaces (bool ReprocessFaces=true)
	Extract element faces. More...

virtual void	ProcessElements ()
	Generate element IDs. More...

virtual void	ProcessComposites ()
	Generate composites. More...

virtual void	ClearElementLinks ()

Static Public Member Functions
static boost::shared_ptr< Module >	create (FieldSharedPtr f)
	Creates an instance of this class. More...

Static Public Member Functions inherited from Nektar::Utilities::ProcessEquiSpacedOutput
static boost::shared_ptr< Module >	create (FieldSharedPtr f)
	Creates an instance of this class. More...

Static Public Attributes
static ModuleKey	className

Static Public Attributes inherited from Nektar::Utilities::ProcessEquiSpacedOutput
static ModuleKey	className

Protected Member Functions
	ProcessIsoContour ()

void	ResetFieldPts (vector< IsoSharedPtr > &iso)

void	SetupIsoFromFieldPts (vector< IsoSharedPtr > &isovec)

Protected Member Functions inherited from Nektar::Utilities::ProcessEquiSpacedOutput
	ProcessEquiSpacedOutput ()

void	SetupEquiSpacedField (void)

void	SetHomogeneousConnectivity (void)

void	GenOrthoModes (int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)

Protected Member Functions inherited from Nektar::Utilities::Module
	Module ()

void	ReorderPrisms (PerMap &perFaces)
	Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...

void	PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)

Private Member Functions
vector< IsoSharedPtr >	ExtractContour (const int fieldid, const NekDouble val)

Additional Inherited Members
Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr	m_f
	Field object. More...

map< string, ConfigOption >	m_config
	List of configuration values. More...

bool	m_requireEquiSpaced

MeshSharedPtr	m_mesh
	Mesh object. More...

std::map< std::string, ConfigOption >	m_config
	List of configuration values. More...

Detailed Description

This processing module extracts an isocontour.

Definition at line 227 of file ProcessIsoContour.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessIsoContour::ProcessIsoContour ( FieldSharedPtr f )

Definition at line 60 of file ProcessIsoContour.cpp.

References Nektar::Utilities::Module::m_config.

                                                      :
     ProcessEquiSpacedOutput(f)
 {
 
     m_config["fieldstr"]           = ConfigOption(false, "NotSet",
                                         "string of isocontour to be extracted");
     m_config["fieldname"]          = ConfigOption(false, "isocon",
                                         "name for isocontour if fieldstr "
                                         "specified, default is isocon");
 
     m_config["fieldid"]            = ConfigOption(false, "NotSet",
                                         "field id to extract");
 
     m_config["fieldvalue"]         = ConfigOption(false, "NotSet",
                                         "field value to extract");
 
     m_config["globalcondense"]     = ConfigOption(true, "NotSet",
                                         "Globally condense contour to unique "
                                         "values");
 
     m_config["smooth"]             = ConfigOption(true, "NotSet",
                                         "Smooth isocontour (might require "
                                                   "globalcondense)");
 
     m_config["smoothiter"]         = ConfigOption(false, "100",
                                         "Number of smoothing cycle, default = "
                                         "100");
 
     m_config["smoothposdiffusion"] = ConfigOption(false,"0.5",
                                         "Postive diffusion coefficient "
                                         "(0 < lambda < 1), default = 0.5");
 
     m_config["smoothnegdiffusion"] = ConfigOption(false,"0.495",
                                         "Negative diffusion coefficient "
                                         "(0 < mu < 1), default = 0.495");
 
     m_config["removesmallcontour"] = ConfigOption(false,"0",
                                         "Remove contours with less than specified number of triangles."
                                          "Only valid with GlobalCondense or Smooth options.");
 
 
 }

Nektar::Utilities::ProcessIsoContour::~ProcessIsoContour ( void )

virtual

Definition at line 103 of file ProcessIsoContour.cpp.

104 {

105 }

Nektar::Utilities::ProcessIsoContour::ProcessIsoContour ( )

inlineprotected

Definition at line 244 of file ProcessIsoContour.h.

244 {};

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessIsoContour::create ( FieldSharedPtr f )

inlinestatic

Creates an instance of this class.

Definition at line 231 of file ProcessIsoContour.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

         {
             return MemoryManager<ProcessIsoContour>::AllocateSharedPtr(f);
         }

vector< IsoSharedPtr > Nektar::Utilities::ProcessIsoContour::ExtractContour	(	const int	fieldid,
		const NekDouble	val
	)

private

Definition at line 355 of file ProcessIsoContour.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::LibUtilities::ePtsTetBlock, Nektar::Utilities::Module::m_f, Nektar::Utilities::ThreeSimilar(), and Nektar::Utilities::TwoPairs().

Referenced by Process().

 {
     vector<IsoSharedPtr> returnval;
 
     int coordim = m_f->m_fieldPts->GetDim();
     int nfields = m_f->m_fieldPts->GetNFields() + coordim;
 
     ASSERTL0(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock,
              "This methods is currently only set up for 3D fields");
     ASSERTL1(coordim + fieldid < nfields,
              "field id is larger than number contained in FieldPts");
     Array<OneD, Array<OneD, NekDouble> > fields;
     m_f->m_fieldPts->GetPts(fields);
 
     Array<OneD, NekDouble> c = fields[coordim + fieldid];
 
     int i, j, k, ii, jj, kk, r, s, n, counter, boolean;
     Array<OneD, Array<OneD, NekDouble> > intfields(nfields);
     intfields[0] = Array<OneD, NekDouble>(5*nfields);
     for(i = 1; i < nfields; ++i)
     {
         intfields[i] = intfields[i-1] + 5;
     }
     Array<OneD, NekDouble> cx = intfields[0];
     Array<OneD, NekDouble> cy = intfields[1];
     Array<OneD, NekDouble> cz = intfields[2];
 
     vector<Array<OneD, int> > ptsConn;
     m_f->m_fieldPts->GetConnectivity(ptsConn);
     for(int zone = 0; zone < ptsConn.size(); ++zone)
     {
         IsoSharedPtr iso;
 
         iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-3);
 
         int nelmt = ptsConn[zone].num_elements()
             /(coordim+1);
 
         Array<OneD, int> conn = ptsConn[zone];
 
         for (n = 0, i = 0; i < nelmt; ++i)
         {
             // check to see if val is between vertex values
             if(!(((c[conn[i*4]]  >val)&&(c[conn[i*4+1]]>val)&&
                   (c[conn[i*4+2]]>val)&&(c[conn[i*4+3]]>val))||
                  ((c[conn[i*4  ]]<val)&&(c[conn[i*4+1]]<val)&&
                   (c[conn[i*4+2]]<val)&&(c[conn[i*4+3]]<val))))
             {
 
                 // loop over all edges and interpolate if
                 // contour is between vertex values
                 for (counter = 0, j=0; j<=2; j++)
                 {
                     for (k=j+1; k<=3; k++)
                     {
                         if (((c[conn[i*4+j]]>=val)&&
                              (val>=c[conn[i*4+k]]))||
                             ((c[conn[i*4+j]]<=val)&&
                              (val<=c[conn[i*4+k]])))
                         {
                             // linear interpolation of fields
                             // (and coords).
                             NekDouble cj = c[conn[i*4+j]];
                             NekDouble ck = c[conn[i*4+k]];
                             NekDouble factor =  (val-cj)/(ck-cj);
 
                             if(fabs(cj-ck) > 1e-12)
                             {
                                 // interpolate coordinates and fields
                                 for(int f = 0; f < nfields; ++f)
                                 {
                                     if(counter == 5)
                                     {
                                         ASSERTL0(false,"Counter is 5");
                                     }
                                     intfields[f][counter] =
                                         fields[f][conn[4*i+j]] +
                                         factor*(fields[f][conn[4*i+k]] -
                                                 fields[f][conn[4*i+j]]);
                                 }
                                 ++counter;
                             }
                         }
                     }
                 }
 
                 switch(counter)
                 {
                 case 3:
                     n+=1;
                     iso->resize_fields(3*n);
 
                     for(j = 0; j < 3; ++j)
                     {
                         iso->set_fields(3*(n-1)+j,intfields,j);
                     }
                     break;
                 case 4:
                     n+=2;
                     iso->resize_fields(3*n);
 
                     for(j = 0; j < 3; ++j)
                     {
                         iso->set_fields(3*(n-2)+j,intfields,j);
                         iso->set_fields(3*(n-1)+j,intfields,j+1);
                     }
                     break;
                 case 5:
                     n+=1;
                     iso->resize_fields(3*n);
 
                     boolean=0;
                     for (ii=0;ii<=2;ii++)
                     {
                         for (jj=ii+1;jj<=3;jj++)
                         {
                             for (kk=jj+1;kk<=4;kk++)
                             {
                                 if((((cx[ii]-cx[jj])==0.0)&&
                                     ((cy[ii]-cy[jj])==0.0)&&
                                     ((cz[ii]-cz[jj])==0.0))&&
                                    (((cx[ii]-cx[kk])==0.0)&&
                                     ((cy[ii]-cy[kk])==0.0)&&
                                     ((cz[ii]-cz[kk])==0.0)))
                                 {
                                     boolean+=1;
                                     ThreeSimilar (ii,jj,kk,r,s);
 
                                     iso->set_fields(3*(n-1)  ,intfields,ii);
                                     iso->set_fields(3*(n-1)+1,intfields,r);
                                     iso->set_fields(3*(n-1)+2,intfields,s);
                                 }
                                 else
                                 {
                                     boolean+=0;
                                 }
                             }
                         }
                     }
 
                     if (boolean==0)
                     {
                         TwoPairs (cx,cy,cz,r);
 
                         iso->set_fields(3*(n-1)  ,intfields,0);
                         iso->set_fields(3*(n-1)+1,intfields,2);
                         iso->set_fields(3*(n-1)+2,intfields,r);
                     }
                     break;
                 }
             }
         }
         iso->set_ntris(n);
 
         // condense the information in this elemental extraction.
         iso->condense();
         returnval.push_back(iso);
     }
 
     return returnval;
 }

void Nektar::Utilities::ProcessIsoContour::Process ( po::variables_map & vm )

virtual

Write mesh to output file.

Reimplemented from Nektar::Utilities::ProcessEquiSpacedOutput.

Definition at line 107 of file ProcessIsoContour.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ExtractContour(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, ResetFieldPts(), Nektar::Utilities::ProcessEquiSpacedOutput::SetupEquiSpacedField(), SetupIsoFromFieldPts(), Nektar::Timer::Start(), Nektar::Timer::Stop(), and Nektar::Timer::TimePerTest().

 {
     Timer timer;
     int rank = m_f->m_comm->GetRank();
     
     if(m_f->m_verbose)
     {
         if(rank == 0)
         {
             cout << "Process Contour extraction..." << endl;
             timer.Start();
         }
     }
 
     vector<IsoSharedPtr> iso;
 
     if(m_f->m_fieldPts.get()) // assume we have read .dat file to directly input dat file.
     {
         SetupIsoFromFieldPts(iso);
     }
     else // extract isocontour from field 
     {
         if(m_f->m_exp.size() == 0)
         {
             return;
         }
 
         // extract all fields to equi-spaced
         SetupEquiSpacedField();
         
         int     fieldid;
         NekDouble value;
         
         if(m_config["fieldstr"].m_beenSet) //generate field of interest
         {
             fieldid = m_f->m_fieldPts->GetNFields();
             
             Array<OneD, NekDouble> pts(m_f->m_fieldPts->GetNpoints());
             
             // evaluate new function
             LibUtilities::AnalyticExpressionEvaluator strEval;
             string varstr = "x y z";
             vector<Array<OneD, const NekDouble> > interpfields;
             
             for(int i = 0; i < m_f->m_fieldPts->GetDim(); ++i)
             {
                 interpfields.push_back(m_f->m_fieldPts->GetPts(i));
             }
             for(int i = 0; i < m_f->m_fieldPts->GetNFields(); ++i)
             {
                 varstr += " " + m_f->m_fieldPts->GetFieldName(i);
                 interpfields.push_back(m_f->m_fieldPts->GetPts(i+3));
             }
             
             int ExprId  = -1;
             std::string str = m_config["fieldstr"].as<string>();
             ExprId = strEval.DefineFunction(varstr.c_str(), str);
             
             strEval.Evaluate(ExprId, interpfields, pts);
             
             // set up field name if provided otherwise called "isocon" from default
             string fieldName = m_config["fieldname"].as<string>();
 
             m_f->m_fieldPts->AddField(pts, fieldName);
         }
         else
         {
             ASSERTL0(m_config["fieldid"].as<string>() != "NotSet", "fieldid must be specified");
             fieldid = m_config["fieldid"].as<int>();
         }
         
         ASSERTL0(m_config["fieldvalue"].as<string>() != "NotSet", "fieldvalue must be specified");
         value   = m_config["fieldvalue"].as<NekDouble>();
 
         iso = ExtractContour(fieldid,value);
 
     }
 
     // Process isocontour
     bool smoothing      = m_config["smooth"].m_beenSet;
     bool globalcondense = m_config["globalcondense"].m_beenSet;
     if(globalcondense)
     {
         int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
         IsoSharedPtr g_iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-3);
 
         g_iso->globalcondense(iso,m_f->m_verbose);
 
 
         iso.clear();
         iso.push_back(g_iso);
     }
 
     if(smoothing)
     {
         int  niter = m_config["smoothiter"].as<int>();
         NekDouble lambda = m_config["smoothposdiffusion"].as<NekDouble>();
         NekDouble mu     = m_config["smoothnegdiffusion"].as<NekDouble>();
         for(int i =0 ; i < iso.size(); ++i)
         {
             iso[i]->smooth(niter,lambda,-mu);
         }
     }
     
     int mincontour = 0; 
     if((mincontour = m_config["removesmallcontour"].as<int>()))
     {
         vector<IsoSharedPtr> new_iso;
 
         if(rank == 0)
         {
             cout << "Identifying separate regions [." << flush ;
         }
         for(int i =0 ; i < iso.size(); ++i)
         {
             iso[i]->separate_regions(new_iso,mincontour,m_f->m_verbose);
         }
        
         if(rank == 0)
         {
             cout << "]" << endl <<  flush ;
         }
 
         // reset iso to new_iso;
         iso = new_iso;
     }
     
     ResetFieldPts(iso);
     
 
     if(m_f->m_verbose)
     {
         if(rank == 0)
         {
             timer.Stop();
             NekDouble cpuTime = timer.TimePerTest(1);
             
             stringstream ss;
             ss << cpuTime << "s";
             cout << "Process Isocontour CPU Time: " << setw(8) << left
                  << ss.str() << endl;
             cpuTime = 0.0;
         }
     }
 }

void Nektar::Utilities::ProcessIsoContour::ResetFieldPts ( vector< IsoSharedPtr > & iso )

protected

Definition at line 520 of file ProcessIsoContour.cpp.

References Nektar::LibUtilities::ePtsTriBlock, Nektar::Utilities::Module::m_f, and npts.

Referenced by Process().

 {
     int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
 
     // set output to triangle block.
     m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
 
     Array<OneD, Array<OneD, NekDouble> > newfields(nfields);
 
     // count up number of points
     int npts = 0;
     for(int i =0; i < iso.size(); ++i)
     {
         npts += iso[i]->get_nvert();
     }
 
     // set up coordinate in new field
     newfields[0] = Array<OneD, NekDouble>(npts);
     newfields[1] = Array<OneD, NekDouble>(npts);
     newfields[2] = Array<OneD, NekDouble>(npts);
 
     int cnt = 0;
     for(int i =0; i < iso.size(); ++i)
     {
         for(int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
         {
             newfields[0][cnt] = iso[i]->get_x(j);
             newfields[1][cnt] = iso[i]->get_y(j);
             newfields[2][cnt] = iso[i]->get_z(j);
         }
     }
 
 
     // set up fields
     for(int f = 0; f < nfields-3; ++f)
     {
         newfields[f+3] = Array<OneD, NekDouble>(npts);
 
         cnt = 0;
         for(int i =0; i < iso.size(); ++i)
         {
             for(int j = 0; j < iso[i]->get_nvert(); ++j,++cnt)
             {
                 newfields[f+3][cnt] = iso[i]->get_fields(f,j);
             }
         }
     }
 
     m_f->m_fieldPts->SetPts(newfields);
 
     // set up connectivity data.
     vector<Array<OneD, int> > ptsConn;
     m_f->m_fieldPts->GetConnectivity(ptsConn);
     cnt = 0;
     ptsConn.clear();
     for(int i =0; i < iso.size(); ++i)
     {
         int ntris = iso[i]->get_ntris();
         Array<OneD, int> conn(ntris*3);
 
         for(int j = 0; j < 3*ntris; ++j)
         {
             conn[j] = cnt + iso[i]->get_vid(j);
         }
         ptsConn.push_back(conn);
         cnt += iso[i]->get_nvert();
     }
     m_f->m_fieldPts->SetConnectivity(ptsConn);
 }

void Nektar::Utilities::ProcessIsoContour::SetupIsoFromFieldPts ( vector< IsoSharedPtr > & isovec )

protected

Definition at line 591 of file ProcessIsoContour.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::ePtsTriBlock, and Nektar::Utilities::Module::m_f.

Referenced by Process().

 {
     ASSERTL0(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock,
              "Assume input is from ePtsTriBlock");
 
     // get information from PtsField
     int dim     = m_f->m_fieldPts->GetDim();
     int nfields = m_f->m_fieldPts->GetNFields() + dim;
     Array<OneD, Array<OneD, NekDouble> > fieldpts;
     m_f->m_fieldPts->GetPts(fieldpts); 
     vector<Array<OneD, int> > ptsConn;
     m_f->m_fieldPts->GetConnectivity(ptsConn);
 
 
     int cnt = 0;
     for(int c = 0; c < ptsConn.size(); ++c)
     {
         // set up single iso with all the information from PtsField
         IsoSharedPtr iso = MemoryManager<Iso>::AllocateSharedPtr(nfields-dim);
         
         int nelmt = 0;
         nelmt = ptsConn[c].num_elements()/3;
 
         iso->set_ntris(nelmt);
         iso->resize_vid(3*nelmt);
         
         // fill in connectivity values. 
         int nvert = 0; 
         for(int i = 0; i < ptsConn[c].num_elements(); ++i)
         {
             int cid = ptsConn[c][i]-cnt;
             iso->set_vid(i,cid);
             nvert = max(cid,nvert);
         }
         nvert++;
         
         iso->set_nvert(nvert);
         iso->resize_fields(nvert);
         
         // fill in points values (including coordinates)
         for(int i = 0; i < nvert; ++i)
         {
             iso->set_fields(i,fieldpts,i+cnt);
         }
         cnt += nvert;
         isovec.push_back(iso);
     }
     
 }

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessIsoContour::className

static

Initial value:

=
    GetModuleFactory().RegisterCreatorFunction(
                        ModuleKey(eProcessModule, "isocontour"),
                        ProcessIsoContour::create,
                        "Extract an isocontour of fieldid variable and at "
                        "value fieldvalue, Optionally fieldstr can be "
                        "specified for a string defiition or smooth for "
                        "smoothing")

Definition at line 235 of file ProcessIsoContour.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation