Nektar::FieldUtils::ProcessBodyFittedVelocity Class Reference

This processing module calculates the wall shear stress and adds it as an extra-field to the output file, and writes it to a surface output file. More...

#include <ProcessBodyFittedVelocity.h>

Inheritance diagram for Nektar::FieldUtils::ProcessBodyFittedVelocity:

Public Member Functions
	ProcessBodyFittedVelocity (FieldSharedPtr f)
virtual	~ProcessBodyFittedVelocity ()
void	GenPntwiseBodyFittedCoordSys (const int targetBndId, const Array< OneD, NekDouble > assistVec, Array< OneD, NekDouble > &distance, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &bfcsDir, const bool isCheckAngle, const NekDouble distTol=1.0e-12, const NekDouble iterTol=1.0e-12, const NekDouble dirTol=1.0e-4, const NekDouble geoTol=1.0e-12)
	At each quadrature point inside the domian, compute the body-fitted coordinate system with respect to the input boundary id. More...
Public Member Functions inherited from Nektar::FieldUtils::ProcessBoundaryExtract
	ProcessBoundaryExtract (FieldSharedPtr f)
virtual	~ProcessBoundaryExtract ()
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 Member Functions inherited from Nektar::FieldUtils::ProcessBoundaryExtract
static std::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::FieldUtils::ProcessBoundaryExtract
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
Protected Member Functions inherited from Nektar::FieldUtils::ProcessBoundaryExtract
virtual void	v_Process (po::variables_map &vm) override
virtual std::string	v_GetModuleName () override
virtual std::string	v_GetModuleDescription () override
virtual ModulePriority	v_GetModulePriority () override
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
NekDouble	PntToBndElmtPntDistance (const Array< OneD, Array< OneD, NekDouble > > &pts, const int pId, const Array< OneD, Array< OneD, NekDouble > > &bndPts)
	Compute the local coordinate for the nearest point on the given 2D boundary element to the input point. More...
bool	LocCoordForNearestPntOnBndElmt_2D (const Array< OneD, const NekDouble > &inGloCoord, SpatialDomains::GeometrySharedPtr bndGeom, const Array< OneD, Array< OneD, NekDouble > > &pts, Array< OneD, NekDouble > &locCoord, Array< OneD, NekDouble > &gloCoord, NekDouble &dist, const NekDouble iterTol=1.0e-12, const int iterMax=51)
	Compute the local coordinate for the nearest point on the given 2D boundary element to the input point. More...
bool	LocCoordForNearestPntOnBndElmt (const Array< OneD, const NekDouble > &inGloCoord, SpatialDomains::GeometrySharedPtr bndGeom, Array< OneD, NekDouble > &locCoord, Array< OneD, NekDouble > &gloCoord, NekDouble &dist, const NekDouble iterTol=1.0e-12, const int iterMax=51)
	Compute the local coordinate for the nearest point on the given boundary element to the input point. The locCoord is the position to set up the body-fitted coordinate. This function works as a driver. More...
void	ScaledCrosssProduct (const Array< OneD, NekDouble > &vec1, const Array< OneD, NekDouble > &vec2, Array< OneD, NekDouble > &vec3)
	Compute the normalized cross product for two 2D vectors. vec3 = vec1 x vec2. More...
void	GetVelAndConvertToCartSys (Array< OneD, Array< OneD, NekDouble > > &vel)
	Get velocity and convert to Cartesian system, if it is still in transformed system. It is copied and modified from from ProcessGrad.cpp. More...

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 calculates the wall shear stress and adds it as an extra-field to the output file, and writes it to a surface output file.

Definition at line 50 of file ProcessBodyFittedVelocity.h.

Constructor & Destructor Documentation

ProcessBodyFittedVelocity()

Nektar::FieldUtils::ProcessBodyFittedVelocity::ProcessBodyFittedVelocity

(

FieldSharedPtr

f

)

Definition at line 64 of file ProcessBodyFittedVelocity.cpp.

    : ProcessBoundaryExtract(f)
{
    f->m_writeBndFld =
        false; // turned on in the father class ProcessBoundaryExtract
 
    // bnd is read from father class ProcessBoundaryExtract
    m_config["assistDir"] =
        ConfigOption(false, "0.0,0.0,1.0",
                     "The normal direction of the assistant plane, where \
                             the first body-fitted coordinate direction is located. \
                             The assistance plane is defined in point normal form. \
                             Default = [0,0,1]");
    m_config["checkAngle"] =
        ConfigOption(true, "0", "The flag for checking the distance vector \
                             perpendicular to the nearest element or not. It \
                             does not need to be turned on if geometry is \
                             smooth and there is no special consideration of \
                             the body-fitted coordinate system. It would be \
                             helpful to set up a desired coordinate system \
                             based on the surface with with sharp corners, \
                             such as at the presence of a step or gap. \
                             Default = false.");
    m_config["distTol"] = ConfigOption(
        false, "1.0e-12", "The distance tolerence to properly find out the \
                             nearest boundary element. It works together with \
                             the CHECKANGLE option to set up the body-fitted \
                             coordinate system near the corner. \
                             Default = 1.0e-12.");
    m_config["bfcOut"] = ConfigOption(
        true, "1", "The flag control whether the body-fitted coordinate \
                             system will be output. Default = true.");
}

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

~ProcessBodyFittedVelocity()

Nektar::FieldUtils::ProcessBodyFittedVelocity::~ProcessBodyFittedVelocity ( )

virtual

Definition at line 98 of file ProcessBodyFittedVelocity.cpp.

{
}

Member Function Documentation

create()

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

inlinestatic

Creates an instance of this class.

Definition at line 54 of file ProcessBodyFittedVelocity.h.

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

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

GenPntwiseBodyFittedCoordSys()

void Nektar::FieldUtils::ProcessBodyFittedVelocity::GenPntwiseBodyFittedCoordSys	(	const int	targetBndId,
		const Array< OneD, NekDouble >	assistVec,
		Array< OneD, NekDouble > &	distance,
		Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &	bfcsDir,
		const bool	isCheckAngle,
		const NekDouble	distTol = 1.0e-12,
		const NekDouble	iterTol = 1.0e-12,
		const NekDouble	dirTol = 1.0e-4,
		const NekDouble	geoTol = 1.0e-12
	)

At each quadrature point inside the domian, compute the body-fitted coordinate system with respect to the input boundary id.

Parameters

targetBndId	Target boundary id.
assistVec	Unit assistant vector, the cross product of inward- pointing wall normalId and wihch gives of the main tangential direction of the body-fitted system.
bfcsDir	Pointwise body-fitted coordinate system.
isPerpendicularCondition	Flag for using perpendicular check or not
distTol	Distance tolerence. Used to find the boundary elements for local coordinate iteration.
iterTol	Iteration tolerence. Used to check iteration convergence.
dirTol	Direction tolerencce. Used to check if the inner product of two unit vectors is cloes enough to 1.0.
geoTol	Geometry tolerence. Used as the relative tolerence for local coord and distance absolute tolerence.

Definition at line 506 of file ProcessBodyFittedVelocity.cpp.

{
 
    const int nFields       = m_f->m_variables.size();
    const int nCoordDim     = m_f->m_exp[0]->GetCoordim(0);
    const int nSpaceDim     = nCoordDim + m_f->m_numHomogeneousDir;
    const int nBndLcoordDim = nCoordDim - 1;
 
    // Get expansion list for boundary
    Array<OneD, MultiRegions::ExpListSharedPtr> BndExp(nFields);
    for (int i = 0; i < nFields; ++i)
    {
        BndExp[i] = m_f->m_exp[i]->UpdateBndCondExpansion(targetBndId);
    }
 
    const int nElmts    = m_f->m_exp[0]->GetNumElmts(); // num of inner element
    const int nBndElmts = BndExp[0]->GetNumElmts();     // num of bnd element
 
    // Use the outward-pointing normal at quardrature pts on bnd elmt as dir ref
    Array<OneD, Array<OneD, NekDouble>> normalQ;
    m_f->m_exp[0]->GetBoundaryNormals(targetBndId, normalQ);
 
    // Vars in the loop
    SpatialDomains::GeometrySharedPtr bndGeom;
    StdRegions::StdExpansionSharedPtr bndXmap;
    Array<OneD, Array<OneD, NekDouble>> bndCoeffs(nCoordDim);
    Array<OneD, Array<OneD, NekDouble>> bndPts(nCoordDim);
 
    int nPts, nBndPts;
    NekDouble dist_tmp, dist_bak;
    bool isConverge;
    vector<int> bndElmtIds;
    int bndElmtId, bndElmtId_bak;
    int phys_offset, bnd_phys_offset;
 
    Array<OneD, NekDouble> inGloCoord(3, 0.0); // inner pnt
    Array<OneD, NekDouble> locCoord(nBndLcoordDim), locCoord_tmp(nBndLcoordDim);
    Array<OneD, NekDouble> gloCoord(3, 0.0), gloCoord_tmp(3, 0.0);
    Array<OneD, NekDouble> normal(3), normalRef(3);
    Array<OneD, NekDouble> tangential1(3),
        tangential2(3); // 1 - main, 2 - minor
    Array<OneD, NekDouble> normalChk(3);
    ProcessWallNormalData wnd(m_f); // wallNormalData object to use its routine
 
    // backup for angle check
    Array<OneD, NekDouble> locCoord_bak(nBndLcoordDim), gloCoord_bak(3, 0.0);
 
    //-------------------------------------------------------------------------
    // Loop all the quadrature points of the field, and find out the closest
    // point on the bnd for each of them. The bnd element contains this bnd
    // point is recorded.
 
    // Loop inner element
    for (int eId = 0; eId < nElmts; ++eId) // element id, nElmts
    {
        // Get inner points
        nPts = m_f->m_exp[0]->GetExp(eId)->GetTotPoints();
        Array<OneD, Array<OneD, NekDouble>> pts(nCoordDim);
        for (int i = 0; i < nCoordDim; ++i)
        {
            pts[i] = Array<OneD, NekDouble>(nPts, 0.0);
        }
 
        if (nCoordDim == 2)
        {
            m_f->m_exp[0]->GetExp(eId)->GetCoords(pts[0], pts[1]);
        }
        else if (nCoordDim == 3)
        {
            m_f->m_exp[0]->GetExp(eId)->GetCoords(pts[0], pts[1], pts[2]);
        }
 
        // Get the offset for the current elmt
        phys_offset = m_f->m_exp[0]->GetPhys_Offset(eId);
 
        // Loop points in the inner element
        for (int pId = 0; pId < nPts; ++pId)
        {
 
            // Compute the distance from the pnt (bnd elmt) to the inner pnt
            // Step 1 - Estimate search
            //  - Find the closest quadrature pnt among all the bnd elmt;
            //  - The bnd elmt which contains this pnt is saved for Step 2.
            //  * If the quadtature pnt is on the edge or corner of the bnd
            //    elmt, there will be more than one bnd elmt which contains
            //    this point.
 
            bndElmtIds.clear(); // clear the bnd elmt id vec for a new inner pnt
 
            // Loop bnd elmt to find the elmt to build body-fitted coord sys
            for (int beId = 0; beId < nBndElmts; ++beId)
            {
                // Get geomery and points
                bndGeom = BndExp[0]->GetExp(beId)->GetGeom();
                bndXmap = bndGeom->GetXmap();
                nBndPts = bndXmap->GetTotPoints();
 
                for (int i = 0; i < nCoordDim; ++i)
                {
                    bndPts[i]    = Array<OneD, NekDouble>(nBndPts);
                    bndCoeffs[i] = bndGeom->GetCoeffs(i); // 0/1/2 for x/y/z
                    bndXmap->BwdTrans(bndCoeffs[i], bndPts[i]);
                }
 
                // Compute the distance (estimate)
                dist_tmp = PntToBndElmtPntDistance(pts, pId, bndPts);
 
                if ((dist_tmp < (distance[phys_offset + pId] - geoTol)) &&
                    (dist_tmp < (distance[phys_offset + pId] - distTol)))
                {
                    // Find a closer quadrature point from a different bnd elmt
                    // Update distance and bnd elmt id vector
                    distance[phys_offset + pId] = dist_tmp;
                    bndElmtIds.clear();
                    bndElmtIds.push_back(beId);
                }
                else if ((dist_tmp < (distance[phys_offset + pId] + geoTol)) ||
                         (dist_tmp < (distance[phys_offset + pId] + distTol)))
                {
                    // The same quadrature point from a different bnd elmt is
                    // found for the other time.
                    // Keep the distance and add the bnd elmt id to the vector
                    bndElmtIds.push_back(beId);
                }
 
            } // end of bnd elmt loop
 
            // Step 2 - Accurate search
            //  - Find the accurate distance and locCoord for the nearest pnt
            //    (to the inner point) on the bnd elmt
            //  - This process will be repeated on all the bnd elmts in the bnd
            //    elmt id vector.
 
            // Generate the coordinate for the inner point
            for (int i = 0; i < nCoordDim; ++i)
            {
                inGloCoord[i] = pts[i][pId];
            }
 
            // reset the dist to make sure the condition for if will be
            // satisfied
            distance[phys_offset + pId] = 9999;
            bndElmtId                   = -1;
            dist_bak                    = 9999;
            bndElmtId_bak               = -1;
 
            // Compute the accurate distance and locCoord
            for (int i = 0; i < bndElmtIds.size(); ++i)
            {
                // Iterate on each of the possible bnd elmt to get the locCoord
                // and the smallest dist
                bndGeom = BndExp[0]->GetExp(bndElmtIds[i])->GetGeom();
 
                isConverge = LocCoordForNearestPntOnBndElmt(
                    inGloCoord, bndGeom, locCoord_tmp, gloCoord_tmp, dist_tmp,
                    iterTol, 51);
 
                if (!isConverge)
                {
                    WARNINGL1(false, "Bisection iteration is not converged!!!");
                }
 
                // Perpendicular check (angle check)
                // When activated (distance > geoTol)
                // if the pnt-to-pnt vector is not parallel with the wall
                // normal, this bnd will be skipped. It is helpful to avoid
                // confusion when there are cornors
                if (isCheckAngle && (dist_tmp > geoTol)) // skip pts on the bnd
                {
                    // Generate the scaled vector
                    Vmath::Vcopy(nCoordDim, gloCoord_tmp, 1, normalChk, 1);
                    Vmath::Neg(nCoordDim, normalChk, 1);
                    Vmath::Vadd(nCoordDim, inGloCoord, 1, normalChk, 1,
                                normalChk, 1);
                    Vmath::Smul(nCoordDim,
                                1.0 / sqrt(Vmath::Dot(nCoordDim, normalChk, 1,
                                                      normalChk, 1)),
                                normalChk, 1, normalChk, 1);
 
                    // Get normal based on the locCoord_tmp
                    wnd.GetNormals(bndGeom, locCoord_tmp, normal);
 
                    // Reverse the direction to make it point into the field
                    bnd_phys_offset = BndExp[0]->GetPhys_Offset(bndElmtIds[i]);
                    for (int j = 0; j < nCoordDim; ++j)
                    {
                        normalRef[j] = normalQ[j][bnd_phys_offset];
                    }
 
                    if (Vmath::Dot(nCoordDim, normal, normalRef) > 0.0)
                    {
                        Vmath::Neg(nCoordDim, normal, 1);
                    }
 
                    // Check if the current bnd elmt convers this inner pnt
                    if (Vmath::Dot(nCoordDim, normalChk, 1, normal, 1) <
                        (1 - dirTol))
                    {
                        // If not covered, save the nearest result as backup
                        // Then skip
                        if (dist_tmp < dist_bak)
                        {
                            bndElmtId_bak = bndElmtIds[i];
                            dist_bak      = dist_tmp;
 
                            for (int j = 0; j < nBndLcoordDim; ++j)
                            {
                                locCoord_bak[j] = locCoord_tmp[j];
                            }
 
                            for (int j = 0; j < 3; ++j)
                            {
                                gloCoord_bak[j] = gloCoord_tmp[j];
                            }
                        }
 
                        continue; // Skip current bnd elmt
                    }
                }
 
                // No violation occurs, update the current result
                // This section must be executed at least once if there the
                // angle check is turned off. It means the following
                // (bndElmtId == -1) can only be caused by angle check.
                if (dist_tmp < distance[phys_offset + pId])
                {
                    bndElmtId                   = bndElmtIds[i];
                    distance[phys_offset + pId] = dist_tmp;
 
                    for (int j = 0; j < nBndLcoordDim; ++j)
                    {
                        locCoord[j] = locCoord_tmp[j];
                    }
 
                    for (int j = 0; j < 3; ++j)
                    {
                        gloCoord[j] = gloCoord_tmp[j];
                    }
                }
            }
 
            // Check if the bnd elmt is found. If not, use nearest one.
            if (bndElmtId == -1)
            {
                if (bndElmtIds.size() == 0)
                {
                    ASSERTL0(false, "No boundary element is found to be the \
                                     possible nearest element. Please check.");
                }
 
                distance[phys_offset + pId] = dist_bak;
                bndElmtId                   = bndElmtId_bak;
 
                for (int i = 0; i < nBndLcoordDim; ++i)
                {
                    locCoord[i] = locCoord_bak[i];
                }
 
                for (int i = 0; i < 3; ++i)
                {
                    gloCoord[i] = gloCoord_bak[i];
                }
 
                WARNINGL1(false,
                          "The boundary element is not found under given \
                                  tolerence. Please check and reset the \
                                  tolerence, or just turn off the angle check. \
                                  If you would like to continue, the nearest \
                                  boundary element is used.");
            }
 
            // Get the wall normal using the function in wallNormalData class
            bndGeom = BndExp[0]->GetExp(bndElmtId)->GetGeom();
            wnd.GetNormals(bndGeom, locCoord, normal);
 
            // Get ref normal for outward-point dir
            bnd_phys_offset = BndExp[0]->GetPhys_Offset(bndElmtId);
            for (int i = 0; i < nCoordDim; ++i)
            {
                normalRef[i] = normalQ[i][bnd_phys_offset];
            }
 
            // Correct the normal direction to make sure it points inward
            if (Vmath::Dot(nCoordDim, normal, normalRef) > 0.0)
            {
                Vmath::Neg(nCoordDim, normal, 1);
            }
 
            // The main tagential direction is defined to be overlapped with the
            // intersection line of the tangantial plane and the assistant
            // plane, whose normal is the input assistDir. Both plane is defined
            // using the point normal form, where the point is the nearest point
            // on the boundary (gloCoord), and the vectors are the normal
            // (pointing out) and the assistant direction (assistDir).
            // Therefore, the main tangantial direction can be obtained by the
            // cross product of the two vectors, since the two vectors has the
            // same starting point and the tangantial direction is perpendicular
            // to both of them.
            ScaledCrosssProduct(normal, assistVec, tangential1);
 
            // Save the direction vectors
            for (int i = 0; i < nSpaceDim; ++i)
            {
                bfcsDir[0][i][phys_offset + pId] = tangential1[i];
                bfcsDir[1][i][phys_offset + pId] = normal[i];
            }
 
            if (nSpaceDim == 3)
            {
                ScaledCrosssProduct(tangential1, normal, tangential2);
                for (int i = 0; i < 3; ++i)
                {
                    bfcsDir[2][i][phys_offset + pId] = tangential2[i];
                }
            }
 
        } // end of inner pts loop
 
    } // end of inner elmt loop
}

References ASSERTL0, Vmath::Dot(), Nektar::FieldUtils::ProcessWallNormalData::GetNormals(), LocCoordForNearestPntOnBndElmt(), Nektar::FieldUtils::Module::m_f, Vmath::Neg(), PntToBndElmtPntDistance(), ScaledCrosssProduct(), Vmath::Smul(), tinysimd::sqrt(), Vmath::Vadd(), Vmath::Vcopy(), and WARNINGL1.

Referenced by v_Process().

GetVelAndConvertToCartSys()

void Nektar::FieldUtils::ProcessBodyFittedVelocity::GetVelAndConvertToCartSys ( Array< OneD, Array< OneD, NekDouble > > &  vel )

private

Get velocity and convert to Cartesian system, if it is still in transformed system. It is copied and modified from from ProcessGrad.cpp.

Definition at line 837 of file ProcessBodyFittedVelocity.cpp.

{
    int spacedim   = m_f->m_exp[0]->GetCoordim(0) + m_f->m_numHomogeneousDir;
    int nfields    = m_f->m_variables.size();
    int npoints    = m_f->m_exp[0]->GetNpoints();
    int var_offset = 0;
 
    // Check type of the fields and set variable offset
    vector<string> vars = m_f->m_exp[0]->GetSession()->GetVariables();
    if (boost::iequals(vars[0], "rho") && boost::iequals(vars[1], "rhou"))
    {
        var_offset = spacedim + 2;
    }
 
    // Get mapping
    GlobalMapping::MappingSharedPtr mapping = ProcessMapping::GetMapping(m_f);
 
    // Get velocity and convert to Cartesian system,
    //      if it is still in transformed system
    if (m_f->m_fieldMetaDataMap.count("MappingCartesianVel"))
    {
        if (m_f->m_fieldMetaDataMap["MappingCartesianVel"] == "False")
        {
            // Initialize arrays and copy velocity
            for (int i = 0; i < spacedim; ++i)
            {
                vel[i] = Array<OneD, NekDouble>(npoints);
                if (m_f->m_exp[0]->GetWaveSpace())
                {
                    m_f->m_exp[0]->HomogeneousBwdTrans(
                        npoints, m_f->m_exp[var_offset + i]->GetPhys(), vel[i]);
                }
                else
                {
                    Vmath::Vcopy(npoints, m_f->m_exp[var_offset + i]->GetPhys(),
                                 1, vel[i], 1);
                }
            }
            // Convert velocity to cartesian system
            mapping->ContravarToCartesian(vel, vel);
            // Convert back to wavespace if necessary
            if (m_f->m_exp[0]->GetWaveSpace())
            {
                for (int i = 0; i < spacedim; ++i)
                {
                    m_f->m_exp[0]->HomogeneousFwdTrans(npoints, vel[i], vel[i]);
                }
            }
        }
        else
        {
            for (int i = 0; i < spacedim; ++i)
            {
                vel[i] = Array<OneD, NekDouble>(npoints);
                Vmath::Vcopy(npoints, m_f->m_exp[var_offset + i]->GetPhys(), 1,
                             vel[i], 1);
            }
        }
    }
    else
    {
        for (int i = 0; i < spacedim && i < nfields; ++i)
        {
            vel[i] = Array<OneD, NekDouble>(npoints);
            Vmath::Vcopy(npoints, m_f->m_exp[var_offset + i]->GetPhys(), 1,
                         vel[i], 1);
        }
    }
}

References Nektar::FieldUtils::ProcessMapping::GetMapping(), Nektar::FieldUtils::Module::m_f, Nektar::GlobalMapping::MappingSharedPtr, and Vmath::Vcopy().

Referenced by v_Process().

LocCoordForNearestPntOnBndElmt()

bool Nektar::FieldUtils::ProcessBodyFittedVelocity::LocCoordForNearestPntOnBndElmt ( const Array< OneD, const NekDouble > &  inGloCoord, SpatialDomains::GeometrySharedPtr  bndGeom, Array< OneD, NekDouble > &  locCoord, Array< OneD, NekDouble > &  gloCoord, NekDouble &  dist, const NekDouble  iterTol = 1.0e-12, const int  iterMax = 51  )

private

Compute the local coordinate for the nearest point on the given boundary element to the input point. The locCoord is the position to set up the body-fitted coordinate. This function works as a driver.

Parameters

inGloCoord	Global coordinate for the input point.
bndGeom	Geometry of the boundary element to search from.
locCoord	Local coordinate of the result.
gloCoord	Global coordinate of the result.
dist	Distance from the input point to the boundary element.
iterTol	Iteration tolerence.
iterMax	Max iteration steps.

Definition at line 430 of file ProcessBodyFittedVelocity.cpp.

{
    bool isConverge     = false;
    const int nCoordDim = m_f->m_exp[0]->GetCoordim(0); // =2 for 2.5D cases
 
    StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
    int nBndPts                               = bndXmap->GetTotPoints();
 
    Array<OneD, Array<OneD, const NekDouble>> bndCoeffs(nCoordDim);
    Array<OneD, Array<OneD, NekDouble>> bndPts(nCoordDim);
    for (int i = 0; i < nCoordDim; ++i)
    {
        bndPts[i]    = Array<OneD, NekDouble>(nBndPts);
        bndCoeffs[i] = bndGeom->GetCoeffs(i);
        bndXmap->BwdTrans(bndCoeffs[i], bndPts[i]);
    }
 
    // Compute the locCoord for the pnt on the bnd elmt
    if (nCoordDim == 2)
    {
        // Bisection search
        isConverge = LocCoordForNearestPntOnBndElmt_2D(
            inGloCoord, bndGeom, bndPts, locCoord, gloCoord, dist, iterTol,
            iterMax);
    }
    else
    {
        // Bi-bisection search
        ASSERTL0(false, "Not available at the moment.");
    }
 
    return isConverge;
}

References ASSERTL0, LocCoordForNearestPntOnBndElmt_2D(), and Nektar::FieldUtils::Module::m_f.

Referenced by GenPntwiseBodyFittedCoordSys().

LocCoordForNearestPntOnBndElmt_2D()

bool Nektar::FieldUtils::ProcessBodyFittedVelocity::LocCoordForNearestPntOnBndElmt_2D ( const Array< OneD, const NekDouble > &  inGloCoord, SpatialDomains::GeometrySharedPtr  bndGeom, const Array< OneD, Array< OneD, NekDouble > > &  pts, Array< OneD, NekDouble > &  locCoord, Array< OneD, NekDouble > &  gloCoord, NekDouble &  dist, const NekDouble  iterTol = 1.0e-12, const int  iterMax = 51  )

private

Compute the local coordinate for the nearest point on the given 2D boundary element to the input point.

Parameters

inGloCoord	Global coordinate for the input point.
bndGeom	Geometry of the boundary element to search from.
pts	Global coordinate of the quadrature points in the boundary element.
locCoord	Local coordinate of the result.
gloCoord	Global coordinate of the result.
dist	Distance from the input point to the boundary element.
iterTol	Iteration tolerence.
iterMax	Max iteration steps.

Definition at line 356 of file ProcessBodyFittedVelocity.cpp.

{
 
    // Initial settings
    Array<OneD, NekDouble> etaLR(2);      // range [-1,1]
    etaLR[0] = -1.0;                      // left
    etaLR[1] = 1.0;                       // right
    NekDouble tmpLx, tmpLy, tmpRx, tmpRy; // tmp values for L/R
    NekDouble distL2, distR2;             // distance square
 
    StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
    locCoord[0]                               = -2.0;
    int cnt                                   = 0;
    bool isConverge                           = false;
    while (cnt < iterMax)
    {
        tmpLx = bndXmap->PhysEvaluate(etaLR, pts[0]);
        tmpLy = bndXmap->PhysEvaluate(etaLR, pts[1]);
        tmpRx = bndXmap->PhysEvaluate(etaLR + 1, pts[0]);
        tmpRy = bndXmap->PhysEvaluate(etaLR + 1, pts[1]);
 
        distL2 = (tmpLx - inGloCoord[0]) * (tmpLx - inGloCoord[0]) +
                 (tmpLy - inGloCoord[1]) * (tmpLy - inGloCoord[1]);
        distR2 = (tmpRx - inGloCoord[0]) * (tmpRx - inGloCoord[0]) +
                 (tmpRy - inGloCoord[1]) * (tmpRy - inGloCoord[1]);
 
        if (distL2 >= distR2)
        {
            etaLR[0] = 0.5 * (etaLR[0] + etaLR[1]);
        }
        else
        {
            etaLR[1] = 0.5 * (etaLR[0] + etaLR[1]);
        }
 
        if ((etaLR[1] - etaLR[0]) < iterTol)
        {
            locCoord[0] = 0.5 * (etaLR[0] + etaLR[1]);
            gloCoord[0] = 0.5 * (tmpLx + tmpRx);
            gloCoord[1] = 0.5 * (tmpLy + tmpRy);
            dist        = sqrt(0.5 * (distL2 + distR2));
            isConverge  = true;
            break;
        }
 
        ++cnt;
    }
 
    // Warning if failed
    if (cnt >= iterMax)
    {
        WARNINGL1(false, "Bisection iteration is not converged");
    }
 
    return isConverge;
}

References tinysimd::sqrt(), and WARNINGL1.

Referenced by LocCoordForNearestPntOnBndElmt().

PntToBndElmtPntDistance()

NekDouble Nektar::FieldUtils::ProcessBodyFittedVelocity::PntToBndElmtPntDistance ( const Array< OneD, Array< OneD, NekDouble > > &  pts, const int  pId, const Array< OneD, Array< OneD, NekDouble > > &  bndPts  )

private

Compute the local coordinate for the nearest point on the given 2D boundary element to the input point.

Compute the point-to-point distance to find the estimated cloest element.

Parameters

inGloCoord	Global coordinate for the input point.
bndGeom	Geometry of the boundary element to search from.
pts	Global coordinate of the quadrature points in the boundary element.
locCoord	Local coordinate of the result.
gloCoord	Global coordinate of the result.
dist	Distance from the input point to the boundary element.
iterTol	Iteration tolerence.
iterMax	Max iteration steps.
pts	Global coordinate of the quadrature points in the inner element.
pId	Id of the inner point of interest in the current loop.
bndPts	Global coordinate of the quadrature points in the boundary element.

Definition at line 315 of file ProcessBodyFittedVelocity.cpp.

{
    const int nCoordDim = pts.size();
 
    NekDouble dist = 9999, dist_tmp;
 
    for (int i = 0; i < bndPts[0].size(); ++i)
    {
        // Compute distance
        dist_tmp = 0.0;
        for (int j = 0; j < nCoordDim; ++j)
        {
            dist_tmp +=
                (bndPts[j][i] - pts[j][pId]) * (bndPts[j][i] - pts[j][pId]);
        }
        dist_tmp = sqrt(dist_tmp);
 
        if (dist_tmp < dist)
        {
            dist = dist_tmp;
        }
    }
 
    return dist;
}

References tinysimd::sqrt().

Referenced by GenPntwiseBodyFittedCoordSys().

ScaledCrosssProduct()

void Nektar::FieldUtils::ProcessBodyFittedVelocity::ScaledCrosssProduct ( const Array< OneD, NekDouble > &  vec1, const Array< OneD, NekDouble > &  vec2, Array< OneD, NekDouble > &  vec3  )

private

Compute the normalized cross product for two 2D vectors. vec3 = vec1 x vec2.

Definition at line 472 of file ProcessBodyFittedVelocity.cpp.

{
    vec3[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
    vec3[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
    vec3[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
 
    NekDouble coef;
    coef =
        1.0 / sqrt(vec3[0] * vec3[0] + vec3[1] * vec3[1] + vec3[2] * vec3[2]);
 
    vec3[0] = vec3[0] * coef;
    vec3[1] = vec3[1] * coef;
    vec3[2] = vec3[2] * coef;
}

References tinysimd::sqrt().

Referenced by GenPntwiseBodyFittedCoordSys().

v_GetModuleDescription()

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

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::ProcessBoundaryExtract.

Definition at line 72 of file ProcessBodyFittedVelocity.h.

    {
        return "Get velocity in the body-fitted coordinate system.";
    }

v_GetModuleName()

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

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::ProcessBoundaryExtract.

Definition at line 67 of file ProcessBodyFittedVelocity.h.

    {
        return "ProcessBodyFittedVelocity";
    }

v_Process()

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

overrideprotectedvirtual

Write mesh to output file.

Reimplemented from Nektar::FieldUtils::ProcessBoundaryExtract.

Definition at line 117 of file ProcessBodyFittedVelocity.cpp.

{
    ProcessBoundaryExtract::v_Process(vm);
 
    // Get dim to store data
    const int nFields    = m_f->m_variables.size();
    const int nCoordDim  = m_f->m_exp[0]->GetCoordim(0);
    const int nSpaceDim  = nCoordDim + m_f->m_numHomogeneousDir;
    const int nAddFields = nSpaceDim + 1; // distanceToWall, u_bfc, v_bfc, w_bfc
 
    // Tolerence setting
    // To include more elmts for further check
    const NekDouble distTol = m_config["distTol"].as<NekDouble>();
    const NekDouble geoTol  = 1.0e-12; // To check if two pnts are too close
    const NekDouble dirTol  = 1.0e-4;  // To check if two unit vecs are parallel
    const NekDouble iterTol = 1.0e-12; // To check if locCoord iter convergence
 
    const bool isCheckAngle = m_config["checkAngle"].as<bool>();
    const bool isOutputBfc  = m_config["bfcOut"].as<bool>();
 
    // Get the assist vector
    vector<NekDouble> assistDir;
    ASSERTL0(ParseUtils::GenerateVector(m_config["assistDir"].as<string>(),
                                        assistDir),
             "Failed to interpret assistance direction");
 
    Array<OneD, NekDouble> assistVec(3); // Save the assist vector
    for (int i = 0; i < 3; ++i)
    {
        assistVec[i] = (assistDir.size() > i) ? assistDir[i] : 0.0;
    }
    if (nCoordDim == 2)
    {
        assistVec[0] = 0.0;
        assistVec[1] = 0.0;
    }
 
    NekDouble norm = sqrt(Vmath::Dot(3, assistVec, 1, assistVec, 1));
    if (norm < geoTol)
    {
        ASSERTL0(false, "Error. The amplitude of assist vector is smaller than \
                         the geometry tolerence 1.0e-12.");
    }
    Vmath::Smul(3, 1.0 / norm, assistVec, 1, assistVec, 1);
 
    // Get the bnd id
    // We only use the first one [0], check ProcessBoundaryExtract for more info
    SpatialDomains::BoundaryConditions bcs(m_f->m_session,
                                           m_f->m_exp[0]->GetGraph());
    const SpatialDomains::BoundaryRegionCollection bregions =
        bcs.GetBoundaryRegions();
    map<int, int> BndRegionMap;
    int cnt = 0;
    for (auto &breg_it : bregions)
    {
        BndRegionMap[breg_it.first] = cnt++;
    }
    int bnd = BndRegionMap[m_f->m_bndRegionsToWrite[0]];
 
    // Generate the distance array and three arrays for directional vectors
    // bfcsDir[i][j][k]
    //   i - dir vec:   0-main tangential, 1-normal, (2-minor tangential)
    //   j - component: 0-x, 1-y, (2-z)
    //   k - point id
    int npoints = m_f->m_exp[0]->GetNpoints();
    Array<OneD, NekDouble> distance(npoints, 9999);
    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> bfcsDir(nSpaceDim);
    for (int i = 0; i < nSpaceDim; ++i)
    {
        bfcsDir[i] = Array<OneD, Array<OneD, NekDouble>>(nSpaceDim);
 
        for (int j = 0; j < nSpaceDim; ++j)
        {
            bfcsDir[i][j] = Array<OneD, NekDouble>(npoints);
        }
    }
 
    // Key function. At each quadrature point inside the domian, compute the
    // body-fitted coordinate system with respect to the input bnd id
    GenPntwiseBodyFittedCoordSys(bnd, assistVec, distance, bfcsDir,
                                 isCheckAngle, distTol, iterTol, dirTol,
                                 geoTol);
 
    // Compute the velocity
    Array<OneD, Array<OneD, NekDouble>> vel_car(nSpaceDim); // Cartesian
    Array<OneD, Array<OneD, NekDouble>> vel_bfc(nSpaceDim); // body-fiitted
    Array<OneD, NekDouble> vel_tmp(npoints);
 
    for (int i = 0; i < nSpaceDim; ++i)
    {
        vel_bfc[i] = Array<OneD, NekDouble>(npoints, 0.0);
    }
 
    // Get the velocity in Cartesian coordinate system
    GetVelAndConvertToCartSys(vel_car);
 
    // Project the velocity into the body-fitted coordinate system
    for (int i = 0; i < nSpaceDim; ++i) // loop for bfc velocity
    {
        for (int j = 0; j < nSpaceDim; ++j)
        {
            Vmath::Vmul(npoints, vel_car[j], 1, bfcsDir[i][j], 1, vel_tmp, 1);
            Vmath::Vadd(npoints, vel_tmp, 1, vel_bfc[i], 1, vel_bfc[i], 1);
        }
    }
 
    // Add var names
    m_f->m_variables.push_back("distanceToWall");
    m_f->m_variables.push_back("u_bfc");
    m_f->m_variables.push_back("v_bfc");
    if (nSpaceDim == 3)
    {
        m_f->m_variables.push_back("w_bfc");
    }
    else if (nSpaceDim == 1)
    {
        ASSERTL0(false, "Velocity in 1D case is already in the body fitted \
                         coordinate. The dimension should be 2 or 3.");
    }
 
    // Add the new fields (distance + u,v,w_bfc)
    // copy distance
    m_f->m_exp.resize(nFields + nAddFields);
 
    m_f->m_exp[nFields] = m_f->AppendExpList(m_f->m_numHomogeneousDir);
    Vmath::Vcopy(npoints, distance, 1, m_f->m_exp[nFields]->UpdatePhys(), 1);
    m_f->m_exp[nFields]->FwdTransLocalElmt(distance,
                                           m_f->m_exp[nFields]->UpdateCoeffs());
 
    // copy vel_bfc
    for (int i = 1; i < nAddFields; ++i)
    {
        m_f->m_exp[nFields + i] = m_f->AppendExpList(m_f->m_numHomogeneousDir);
        Vmath::Vcopy(npoints, vel_bfc[i - 1], 1,
                     m_f->m_exp[nFields + i]->UpdatePhys(), 1);
        m_f->m_exp[nFields + i]->FwdTransLocalElmt(
            vel_bfc[i - 1], m_f->m_exp[nFields + i]->UpdateCoeffs());
    }
 
    // Output the the body-fitted coordinate into the file
    // The filter will read the coordiante through file
    // bfcsDir[i][j][k]
    //   i - dir vec:   0-main tangential, 1-normal, (2-minor tangential)
    //   j - component: 0-x, 1-y, (2-z)
    //   k - point id
    if (isOutputBfc)
    {
        m_f->m_variables.push_back("bfc_dir0_x");
        m_f->m_variables.push_back("bfc_dir0_y");
        if (nSpaceDim == 3)
        {
            m_f->m_variables.push_back("bfc_dir0_z");
        }
        m_f->m_variables.push_back("bfc_dir1_x");
        m_f->m_variables.push_back("bfc_dir1_y");
        if (nSpaceDim == 3)
        {
            m_f->m_variables.push_back("bfc_dir1_z");
 
            m_f->m_variables.push_back("bfc_dir2_x");
            m_f->m_variables.push_back("bfc_dir2_y");
            m_f->m_variables.push_back("bfc_dir2_z");
        }
 
        int nAddFields_coord = (nSpaceDim == 3) ? 9 : 4;
        m_f->m_exp.resize(nFields + nAddFields + nAddFields_coord);
 
        cnt = 0;
        for (int i = 0; i < nSpaceDim; ++i)
        {
            for (int j = 0; j < nSpaceDim; ++j)
            {
 
                m_f->m_exp[nFields + nAddFields + cnt] =
                    m_f->AppendExpList(m_f->m_numHomogeneousDir);
 
                Vmath::Vcopy(
                    npoints, bfcsDir[i][j], 1,
                    m_f->m_exp[nFields + nAddFields + cnt]->UpdatePhys(), 1);
                m_f->m_exp[nFields + nAddFields + cnt]->FwdTransLocalElmt(
                    bfcsDir[i][j],
                    m_f->m_exp[nFields + nAddFields + cnt]->UpdateCoeffs());
 
                ++cnt;
            }
        }
    }
}

References ASSERTL0, Vmath::Dot(), Nektar::ParseUtils::GenerateVector(), GenPntwiseBodyFittedCoordSys(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), GetVelAndConvertToCartSys(), Nektar::FieldUtils::Module::m_config, Nektar::FieldUtils::Module::m_f, Vmath::Smul(), tinysimd::sqrt(), Nektar::FieldUtils::ProcessBoundaryExtract::v_Process(), Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vmul().

Member Data Documentation

className

ModuleKey Nektar::FieldUtils::ProcessBodyFittedVelocity::className

static

Initial value:

=
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "bodyFittedVelocity"),
        ProcessBodyFittedVelocity::create,
        "Convert the velocity components from the Cartesian coordinate into "
        "the body-fitted coordinate.")

Definition at line 58 of file ProcessBodyFittedVelocity.h.