Nektar::FieldUtils::ProcessIsoContour Class Reference

This processing module extracts an isocontour. More...

#include <ProcessIsoContour.h>

Inheritance diagram for Nektar::FieldUtils::ProcessIsoContour:

Public Member Functions
	ProcessIsoContour (FieldSharedPtr f)
virtual	~ProcessIsoContour ()
Public Member Functions inherited from Nektar::FieldUtils::ProcessModule
	ProcessModule ()
	ProcessModule (FieldSharedPtr p_f)
Public Member Functions inherited from Nektar::FieldUtils::Module
FIELD_UTILS_EXPORT	Module (FieldSharedPtr p_f)
virtual	~Module ()=default
void	Process (po::variables_map &vm)
std::string	GetModuleName ()
std::string	GetModuleDescription ()
const ConfigOption &	GetConfigOption (const std::string &key) const
ModulePriority	GetModulePriority ()
FIELD_UTILS_EXPORT void	RegisterConfig (std::string key, std::string value="")
	Register a configuration option with a module. More...
FIELD_UTILS_EXPORT void	PrintConfig ()
	Print out all configuration options for a module. More...
FIELD_UTILS_EXPORT void	SetDefaults ()
	Sets default configuration options for those which have not been set. More...
FIELD_UTILS_EXPORT void	AddFile (std::string fileType, std::string fileName)
FIELD_UTILS_EXPORT void	EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)

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

Static Public Attributes
static ModuleKey	className

Protected Member Functions
virtual void	v_Process (po::variables_map &vm) override
	Write mesh to output file. More...
virtual std::string	v_GetModuleName () override
virtual std::string	v_GetModuleDescription () override
virtual ModulePriority	v_GetModulePriority () override
	ProcessIsoContour ()
void	ResetFieldPts (std::vector< IsoSharedPtr > &iso)
void	SetupIsoFromFieldPts (std::vector< IsoSharedPtr > &isovec)
Protected Member Functions inherited from Nektar::FieldUtils::Module
	Module ()
virtual void	v_Process (po::variables_map &vm)
virtual std::string	v_GetModuleName ()
virtual std::string	v_GetModuleDescription ()
virtual ModulePriority	v_GetModulePriority ()

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

Additional Inherited Members
Public Attributes inherited from Nektar::FieldUtils::Module
FieldSharedPtr	m_f
	Field object. More...
Protected Attributes inherited from Nektar::FieldUtils::Module
std::map< std::string, ConfigOption >	m_config
	List of configuration values. More...
std::set< std::string >	m_allowedFiles
	List of allowed file formats. More...

Detailed Description

This processing module extracts an isocontour.

Definition at line 225 of file ProcessIsoContour.h.

Constructor & Destructor Documentation

ProcessIsoContour() [1/2]

Nektar::FieldUtils::ProcessIsoContour::ProcessIsoContour

(

FieldSharedPtr

f

)

Definition at line 68 of file ProcessIsoContour.cpp.

                                                     : ProcessModule(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, "0",
                     "Globally condense contour to unique "
                     "values");
 
    m_config["smooth"] = ConfigOption(true, "0",
                                      "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.505",
                     "Negative diffusion coefficient "
                     "(0 < mu < 1), default = 0.505");
 
    m_config["removesmallcontour"] = ConfigOption(
        false, "0",
        "Remove contours with less than specified number of triangles."
        "Only valid with GlobalCondense or Smooth options.");
}

References Nektar::FieldUtils::Module::m_config.

~ProcessIsoContour()

Nektar::FieldUtils::ProcessIsoContour::~ProcessIsoContour ( void  )

virtual

Definition at line 112 of file ProcessIsoContour.cpp.

{
}

ProcessIsoContour() [2/2]

Nektar::FieldUtils::ProcessIsoContour::ProcessIsoContour ( )

inlineprotected

Definition at line 258 of file ProcessIsoContour.h.

{};

Member Function Documentation

create()

static std::shared_ptr< Module > Nektar::FieldUtils::ProcessIsoContour::create ( FieldSharedPtr  f )

inlinestatic

Creates an instance of this class.

Definition at line 229 of file ProcessIsoContour.h.

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

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

ExtractContour()

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

private

Definition at line 373 of file ProcessIsoContour.cpp.

{
    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].size() / (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->ResizeFields(3 * n);
 
                        for (j = 0; j < 3; ++j)
                        {
                            iso->SetFields(3 * (n - 1) + j, intfields, j);
                        }
                        break;
                    case 4:
                        n += 2;
                        iso->ResizeFields(3 * n);
 
                        for (j = 0; j < 3; ++j)
                        {
                            iso->SetFields(3 * (n - 2) + j, intfields, j);
                            iso->SetFields(3 * (n - 1) + j, intfields, j + 1);
                        }
                        break;
                    case 5:
                        n += 1;
                        iso->ResizeFields(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->SetFields(3 * (n - 1), intfields,
                                                       ii);
                                        iso->SetFields(3 * (n - 1) + 1,
                                                       intfields, r);
                                        iso->SetFields(3 * (n - 1) + 2,
                                                       intfields, s);
                                    }
                                    else
                                    {
                                        boolean += 0;
                                    }
                                }
                            }
                        }
 
                        if (boolean == 0)
                        {
                            TwoPairs(cx, cy, cz, r);
 
                            iso->SetFields(3 * (n - 1), intfields, 0);
                            iso->SetFields(3 * (n - 1) + 1, intfields, 2);
                            iso->SetFields(3 * (n - 1) + 2, intfields, r);
                        }
                        break;
                }
            }
        }
 
        if (n)
        {
            iso->SetNTris(n);
 
            // condense the information in this elemental extraction.
            iso->Condense();
 
            returnval.push_back(iso);
        }
    }
 
    return returnval;
}

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

Referenced by v_Process().

ResetFieldPts()

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

protected

Definition at line 545 of file ProcessIsoContour.cpp.

{
    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]->GetNVert();
    }
 
    // 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]->GetNVert(); ++j, ++cnt)
        {
            newfields[0][cnt] = iso[i]->GetX(j);
            newfields[1][cnt] = iso[i]->GetY(j);
            newfields[2][cnt] = iso[i]->GetZ(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]->GetNVert(); ++j, ++cnt)
            {
                newfields[f + 3][cnt] = iso[i]->GetFields(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]->GetNTris();
        Array<OneD, int> conn(ntris * 3);
 
        for (int j = 0; j < 3 * ntris; ++j)
        {
            conn[j] = cnt + iso[i]->GetVId(j);
        }
        ptsConn.push_back(conn);
        cnt += iso[i]->GetNVert();
    }
    m_f->m_fieldPts->SetConnectivity(ptsConn);
}

References Nektar::LibUtilities::ePtsTriBlock, and Nektar::FieldUtils::Module::m_f.

Referenced by v_Process().

SetupIsoFromFieldPts()

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

protected

Definition at line 615 of file ProcessIsoContour.cpp.

{
    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].size() / 3;
 
        iso->SetNTris(nelmt);
        iso->ResizeVId(3 * nelmt);
 
        // fill in connectivity values.
        int nvert = 0;
        for (int i = 0; i < ptsConn[c].size(); ++i)
        {
            int cid = ptsConn[c][i] - cnt;
            iso->SetVId(i, cid);
            nvert = max(cid, nvert);
        }
        nvert++;
 
        iso->SetNVert(nvert);
        iso->ResizeFields(nvert);
 
        // fill in points values (including coordinates)
        for (int i = 0; i < nvert; ++i)
        {
            iso->SetFields(i, fieldpts, i + cnt);
        }
        cnt += nvert;
        isovec.push_back(iso);
    }
}

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

Referenced by v_Process().

v_GetModuleDescription()

virtual std::string Nektar::FieldUtils::ProcessIsoContour::v_GetModuleDescription ( )

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 247 of file ProcessIsoContour.h.

    {
        return "Extracting contour";
    }

v_GetModuleName()

virtual std::string Nektar::FieldUtils::ProcessIsoContour::v_GetModuleName ( )

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 242 of file ProcessIsoContour.h.

    {
        return "ProcessIsoContour";
    }

v_GetModulePriority()

virtual ModulePriority Nektar::FieldUtils::ProcessIsoContour::v_GetModulePriority ( )

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 252 of file ProcessIsoContour.h.

    {
        return eModifyPts;
    }

References Nektar::FieldUtils::eModifyPts.

v_Process()

void Nektar::FieldUtils::ProcessIsoContour::v_Process ( po::variables_map &  vm )

overrideprotectedvirtual

Write mesh to output file.

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 116 of file ProcessIsoContour.cpp.

{
    m_f->SetUpExp(vm);
 
    bool verbose = (m_f->m_verbose && m_f->m_comm->TreatAsRankZero());
 
    vector<IsoSharedPtr> iso;
 
    ASSERTL0(m_f->m_fieldPts.get(), "Should have m_fieldPts for IsoContour.");
 
    if (m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock)
    {
        // assume we have read .dat file to directly input dat file.
        if (verbose)
        {
            cout << "\t Process read iso from Field Pts" << endl;
        }
 
        SetupIsoFromFieldPts(iso);
    }
    else if (m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock)
    {
        if (m_config["fieldstr"].m_beenSet)
        {
            string fieldName = m_config["fieldname"].as<string>();
            m_f->m_variables.push_back(fieldName);
        }
 
        if (m_f->m_fieldPts->GetNpoints() == 0)
        {
            return;
        }
 
        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::Interpreter 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);
    }
    else
    {
        ASSERTL0(false, "PtsType not supported for isocontour.");
    }
 
    // Process isocontour
    bool smoothing      = m_config["smooth"].as<bool>();
    bool globalcondense = m_config["globalcondense"].as<bool>();
    if (globalcondense)
    {
        if (verbose)
        {
            cout << "\t Process global condense ..." << endl;
        }
        int nfields = m_f->m_fieldPts->GetNFields() + m_f->m_fieldPts->GetDim();
        IsoSharedPtr g_iso = MemoryManager<Iso>::AllocateSharedPtr(nfields - 3);
 
        g_iso->GlobalCondense(iso, verbose);
 
        iso.clear();
        iso.push_back(g_iso);
    }
 
    if (smoothing)
    {
        LibUtilities::Timer timersm;
 
        if (verbose)
        {
            cout << "\t Process Contour smoothing ..." << endl;
            timersm.Start();
        }
 
        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);
        }
 
        if (verbose)
        {
            timersm.Stop();
            NekDouble cpuTime = timersm.TimePerTest(1);
 
            stringstream ss;
            ss << cpuTime << "s";
            cout << "\t Process smooth CPU Time: " << setw(8) << left
                 << ss.str() << endl;
            cpuTime = 0.0;
        }
    }
 
    int mincontour = 0;
    if ((mincontour = m_config["removesmallcontour"].as<int>()))
    {
        vector<IsoSharedPtr> new_iso;
 
        if (verbose)
        {
            cout << "\t Identifying separate regions [." << flush;
        }
        for (int i = 0; i < iso.size(); ++i)
        {
            iso[i]->SeparateRegions(new_iso, mincontour, m_f->m_verbose);
        }
 
        if (verbose)
        {
            cout << "]" << endl << flush;
        }
 
        // reset iso to new_iso;
        iso = new_iso;
    }
 
    ResetFieldPts(iso);
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::Interpreter::DefineFunction(), Nektar::LibUtilities::ePtsTetBlock, Nektar::LibUtilities::ePtsTriBlock, Nektar::LibUtilities::Interpreter::Evaluate(), ExtractContour(), Nektar::FieldUtils::Module::m_config, Nektar::FieldUtils::Module::m_f, ResetFieldPts(), SetupIsoFromFieldPts(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Nektar::LibUtilities::Timer::TimePerTest().

Member Data Documentation

className

ModuleKey Nektar::FieldUtils::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 233 of file ProcessIsoContour.h.