Nektar::Utilities::ProcessMultiShear Class Reference

This processing module calculates the shear stress metrics and writes it to a surface output file. More...

#include <ProcessMultiShear.h>

Inheritance diagram for Nektar::Utilities::ProcessMultiShear:

Collaboration diagram for Nektar::Utilities::ProcessMultiShear:

Public Member Functions
	ProcessMultiShear (FieldSharedPtr f)
virtual	~ProcessMultiShear ()
virtual void	Process (po::variables_map &vm)
	Write mesh to output file. More...
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 Attributes
static ModuleKey	className

Additional Inherited Members
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)
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 calculates the shear stress metrics and writes it to a surface output file.

Definition at line 50 of file ProcessMultiShear.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessMultiShear::ProcessMultiShear

(

FieldSharedPtr

f

)

Definition at line 57 of file ProcessMultiShear.cpp.

References ASSERTL0, Nektar::Utilities::Module::m_config, and Nektar::Utilities::Module::m_f.

                                                     : ProcessModule(f)
{
    m_config["N"] = ConfigOption(false,"1","Number of chk or fld files");
    m_config["fromfld"] = ConfigOption(false, "NotSet",
                                       "First fld file. First underscore flags position of id in name.");

    ASSERTL0(m_config["fromfld"].as<string>().compare("NotSet") != 0,
             "Need to specify fromfld=file.fld ");

    m_f->m_fldToBnd = false;
}

Nektar::Utilities::ProcessMultiShear::~ProcessMultiShear ( )

virtual

Definition at line 69 of file ProcessMultiShear.cpp.

{
}

Member Function Documentation

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

inlinestatic

Creates an instance of this class.

Definition at line 54 of file ProcessMultiShear.h.

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

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

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

virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 73 of file ProcessMultiShear.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, Nektar::LibUtilities::NullFieldMetaDataMap, Vmath::Smul(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vvtvm(), Vmath::Vvtvp(), and Vmath::Zero().

{
    if (m_f->m_verbose)
    {
        cout << "ProcessMultiShear: Calculating shear stress metrics..." << endl;
    }

    int nstart, i, j, nfields;
    bool wssg = false;
    NekDouble nfld = m_config["N"].as<NekDouble>();
    string fromfld, basename, endname, nstartStr;
    stringstream filename;
    vector<string> infiles(nfld);
    vector< boost::shared_ptr<Field> > m_fromField(nfld);

    // Set up list of input fld files.
    fromfld = m_config["fromfld"].as<string>();
    basename = fromfld.substr(0, fromfld.find_first_of("_")+1);
    filename << fromfld.substr(fromfld.find_first_of("_")+1, fromfld.size());
    filename >> nstart;
    filename.str("");
    filename << nstart;
    filename >> nstartStr;
    filename.str("");
    endname = fromfld.substr(fromfld.find(nstartStr)+nstartStr.size(), fromfld.size());

    for (i=0; i<nfld; ++i)
    {
        stringstream filename;
        filename << basename << i+nstart << endname;
        filename >> infiles[i];
        cout << infiles[i]<<endl;
    }

    for ( i = 0; i<nfld; ++i)
    {
        m_fromField[i] = boost::shared_ptr<Field>(new Field());
        m_fromField[i]->m_session = m_f->m_session;
        m_fromField[i]->m_graph = m_f->m_graph;
        m_fromField[i]->m_fld = MemoryManager<LibUtilities::FieldIO>
            ::AllocateSharedPtr(m_fromField[0]->m_session->GetComm());
    }

    //Import all fld files.
    for (i = 0; i < nfld; ++i)
    {
        if(m_f->m_exp.size())
        {
            // Set up ElementGIDs in case of parallel processing
            Array<OneD,int> ElementGIDs(m_f->m_exp[0]->GetExpSize());
            for (j = 0; j < m_f->m_exp[0]->GetExpSize(); ++j)
            {
                ElementGIDs[j] = m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
            }
            m_fromField[i]->m_fld->Import(infiles[i],m_fromField[i]->m_fielddef,
                                          m_fromField[i]->m_data,
                                          LibUtilities::NullFieldMetaDataMap,
                                          ElementGIDs);
        }
        else
        {
            m_fromField[i]->m_fld->Import(infiles[i],m_fromField[i]->m_fielddef,
                                          m_fromField[i]->m_data,
                                          LibUtilities::NullFieldMetaDataMap);
        }

        nfields    = m_fromField[i]->m_fielddef[0]->m_fields.size();
        int NumHomogeneousDir = m_fromField[i]->m_fielddef[0]->m_numHomogeneousDir;

        if (nfields == 5 || nfields == 7)
        {
            wssg = true;
        }

        // Set up Expansion information to use mode order from field
        m_fromField[i]->m_graph->SetExpansions(m_fromField[i]->m_fielddef);

        //Set up expansions, and extract data.
        m_fromField[i]->m_exp.resize(nfields);
        m_fromField[i]->m_exp[0] = m_fromField[i]->SetUpFirstExpList(NumHomogeneousDir,true);

        for(j = 1; j < nfields; ++j)
        {
            m_fromField[i]->m_exp[j] = m_f->AppendExpList(NumHomogeneousDir);
        }

        for (j = 0; j < nfields; ++j)
        {
            for (int k = 0; k < m_fromField[i]->m_data.size(); ++k)
            {
                m_fromField[i]->m_exp[j]->ExtractDataToCoeffs(
                    m_fromField[i]->m_fielddef[k],
                    m_fromField[i]->m_data[k],
                    m_fromField[i]->m_fielddef[k]->m_fields[j],
                    m_fromField[i]->m_exp[j]->UpdateCoeffs());
            }
            m_fromField[i]->m_exp[j]->BwdTrans(m_fromField[i]->m_exp[j]->GetCoeffs(),
                                               m_fromField[i]->m_exp[j]->UpdatePhys());
        }
    }


    int spacedim   = m_f->m_graph->GetSpaceDimension();
    if ((m_fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
        (m_fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 2)
    {
        spacedim = 3;
    }


    int nout = 5; // TAWSS, OSI, transWSS, AFI, CFI (optional WSSG)
    if (wssg) { nout = 6; }

    int npoints = m_fromField[0]->m_exp[0]->GetNpoints();
    Array<OneD, Array<OneD, NekDouble> > normTemporalMeanVec(spacedim),normCrossDir(spacedim),  outfield(nout), dtemp(spacedim);
    Array<OneD, NekDouble> TemporalMeanMag(npoints,0.0), DotProduct(npoints,0.0), temp(npoints,0.0);


    for (i = 0; i < spacedim; ++i)
    {
        normTemporalMeanVec[i] = Array<OneD, NekDouble>(npoints);
        normCrossDir[i] = Array<OneD, NekDouble>(npoints);
        dtemp[i] = Array<OneD, NekDouble>(npoints);
        Vmath::Zero(npoints, normTemporalMeanVec[i],1);
        Vmath::Zero(npoints, normCrossDir[i],1);
    }
    
    for (i = 0; i < nout; ++i)
    {
        outfield[i] = Array<OneD, NekDouble>(npoints, 0.0);
    }

    // -----------------------------------------------------
    // Compute temporal average wall shear stress vector,
    // it's spatial average, and normalise it.
    for (i = 0; i < nfld; ++i)
    {
        for (j = 0; j < spacedim; ++j)
        {
            Vmath::Vadd(npoints,m_fromField[i]->m_exp[j]->GetPhys(), 1, normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1); 
        }
    }

    for (i = 0; i < spacedim; ++i)
    {
        Vmath::Smul(npoints, 1.0/nfld, normTemporalMeanVec[i], 1, normTemporalMeanVec[i], 1);
        Vmath::Vvtvp(npoints, normTemporalMeanVec[i], 1, normTemporalMeanVec[i], 1,
                     TemporalMeanMag, 1, TemporalMeanMag, 1);
    }

    Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);

    for (i = 0; i < spacedim; ++i)
    {
        Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1, normTemporalMeanVec[i], 1);
    }
    //---------------------------------------------------

    if (wssg) //cross product with normals to obtain direction parallel to temporal mean vector. 
    {
        Vmath::Vmul(npoints,m_fromField[0]->m_exp[nfields-1]->GetPhys(),1,normTemporalMeanVec[1],1,normCrossDir[0],1);
        Vmath::Vvtvm(npoints,m_fromField[0]->m_exp[nfields-2]->GetPhys(),1,normTemporalMeanVec[2],1,normCrossDir[0],1,
                     normCrossDir[0],1);

        Vmath::Vmul(npoints,m_fromField[0]->m_exp[nfields-3]->GetPhys(),1,normTemporalMeanVec[2],1,normCrossDir[1],1);
        Vmath::Vvtvm(npoints,m_fromField[0]->m_exp[nfields-1]->GetPhys(),1,normTemporalMeanVec[0],1,normCrossDir[1],1,
                     normCrossDir[1],1);
        
        Vmath::Vmul(npoints,m_fromField[0]->m_exp[nfields-2]->GetPhys(),1,normTemporalMeanVec[0],1,normCrossDir[2],1);
        Vmath::Vvtvm(npoints,m_fromField[0]->m_exp[nfields-3]->GetPhys(),1,normTemporalMeanVec[1],1,normCrossDir[2],1,
                     normCrossDir[2],1);
    }


    // Compute tawss, trs,  osi, taafi, tacfi, WSSG.
    for (i = 0; i < nfld; ++i)
    {
        for (j = 0; j < spacedim; ++j)
        {
            Vmath::Vvtvp(npoints, m_fromField[i]->m_exp[j]->GetPhys(), 1, normTemporalMeanVec[j], 1,
                         DotProduct, 1, DotProduct, 1);
        }
        
        //TAWSS
        Vmath::Vadd(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, outfield[0], 1, outfield[0], 1);
        
        //transWSS
        Vmath::Vmul(npoints, DotProduct, 1, DotProduct, 1, temp,1);
        Vmath::Vvtvm(npoints, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, 
                     m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1, temp, 1);
        
        for (j = 0; j < npoints; ++j)
        {
            if(temp[j] > 0.0)
            {
                outfield[1][j] = outfield[1][j] + sqrt(temp[j]);
            }
        }
        
        //TAAFI
        Vmath::Vdiv(npoints, DotProduct, 1, m_fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1);
        Vmath::Vadd(npoints, temp, 1, outfield[3], 1, outfield[3], 1);
        
        //TACFI
        for (j = 0; j < npoints; ++j)
        {
            temp[j] = 1 - temp[j]*temp[j];
            if(temp[j] > 0.0)
            {
                outfield[4][j] = outfield[4][j] + sqrt(temp[j]);
            }
        }
        
        // WSSG
        if (wssg)
        {
            //parallel component:
            Vmath::Zero(npoints, temp,1);

            m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct,dtemp[0],dtemp[1],dtemp[2]);
            for(j=0; j<spacedim; j++)
            {
                Vmath::Vvtvp(npoints,dtemp[j],1,normTemporalMeanVec[j],1,temp,1,temp,1); 
            }
            Vmath::Vmul(npoints,temp,1,temp,1,temp,1);


            //perpendicular component:
            Vmath::Zero(npoints, DotProduct,1);
           
            for(j = 0; j < spacedim; ++j)
            {
                Vmath::Vvtvp(npoints,m_fromField[i]->m_exp[j]->GetPhys(),1,normCrossDir[j],1,
                             DotProduct, 1, DotProduct, 1);
                
            }
            m_fromField[i]->m_exp[0]->PhysDeriv(DotProduct,dtemp[0],dtemp[1],dtemp[2]);
            Vmath::Zero(npoints, DotProduct,1);

            for(j=0; j<spacedim; j++)
            {
                Vmath::Vvtvp(npoints,dtemp[j],1,normCrossDir[j],1,DotProduct,1,DotProduct,1); 
            }
            Vmath::Vvtvp(npoints,DotProduct,1,DotProduct,1,temp,1, temp,1);

            //Final wssg computation
            Vmath::Vsqrt(npoints,temp,1,temp,1);
            Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
        }
        
        
        Vmath::Zero(npoints, DotProduct,1);
    }

    //Divide by nfld
    Vmath::Smul(npoints, 1.0/nfld, outfield[0], 1, outfield[0], 1);
    Vmath::Smul(npoints, 1.0/nfld, outfield[1], 1, outfield[1], 1);
    Vmath::Smul(npoints, 1.0/nfld, outfield[3], 1, outfield[3], 1);
    Vmath::Smul(npoints, 1.0/nfld, outfield[4], 1, outfield[4], 1);

    //OSI
    for (i = 0; i < npoints; ++i)
    {
        outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i]/outfield[0][i]);
    }

    /* TAWSS = sum(wss)/nfld
     * transWSS = sum( sqrt( wss^2 - (wss . normTempMean)^2) )/nfld.
     * OSI = 0.5*(1-TemporalMeanMag/TAWSS)
     * TAAFI = sum(cos)/nfld
     * TACFI = sum(sin)/nfld = sum( sqrt(1-cos^2) )/nfld.
     */

    m_f->m_exp.resize(nout);
    m_f->m_fielddef = m_fromField[0]->m_fielddef;
    m_f->m_exp[0] = m_f->SetUpFirstExpList(m_f->m_fielddef[0]->m_numHomogeneousDir,true);

    for(i = 1; i < nout; ++i)
    {
        m_f->m_exp[i] = m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir);
    }

    m_f->m_fielddef[0]->m_fields.resize(nout);
    m_f->m_fielddef[0]->m_fields[0] = "TAWSS";
    m_f->m_fielddef[0]->m_fields[1] = "transWSS";
    m_f->m_fielddef[0]->m_fields[2] = "OSI";
    m_f->m_fielddef[0]->m_fields[3] = "TAAFI";
    m_f->m_fielddef[0]->m_fields[4] = "TACFI";

    if (wssg)
    { 
        m_f->m_fielddef[0]->m_fields[5] = "|WSSG|"; 
        Vmath::Smul(npoints, 1.0/nfld, outfield[5], 1, outfield[5], 1);
    }


    for(i = 0; i < nout; ++i)
    {
        m_f->m_exp[i]->FwdTrans(outfield[i],
                                m_f->m_exp[i]->UpdateCoeffs());
        m_f->m_exp[i]->BwdTrans(m_f->m_exp[i]->GetCoeffs(),
                                m_f->m_exp[i]->UpdatePhys());
    }


    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
        = m_fromField[0]->m_exp[0]->GetFieldDefinitions();
    std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());

    for( i = 0; i < nout; ++i)
    {
        for ( j = 0; j < FieldDef.size(); ++j)
        {
            FieldDef[j]->m_fields.push_back(m_f->m_fielddef[0]->m_fields[i]);
            m_f->m_exp[i]->AppendFieldData(FieldDef[j], FieldData[j]);
        }
    }

    m_f->m_fielddef = FieldDef;
    m_f->m_data     = FieldData;
}

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessMultiShear::className

static

Initial value:

=
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "shear"),
        ProcessMultiShear::create, "Computes shear stress metrics.")

Definition at line 57 of file ProcessMultiShear.h.