Nektar::Utilities::ProcessIsoContour Class Reference

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 (string key, 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, set< int > &prismsDone, 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...

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.

{
}

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

inlineprotected

Definition at line 244 of file ProcessIsoContour.h.

{};

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.