Nektar::SpatialDomains::MeshGraphXml Class Reference

#include <MeshGraphXml.h>

Inheritance diagram for Nektar::SpatialDomains::MeshGraphXml:

Public Member Functions
	MeshGraphXml ()
virtual	~MeshGraphXml ()
void	WriteXMLGeometry (std::string outname, std::vector< std::set< unsigned int >> elements, std::vector< unsigned int > partitions)
Public Member Functions inherited from Nektar::SpatialDomains::MeshGraph
	MeshGraph ()
virtual	~MeshGraph ()
void	WriteGeometry (std::string &outfilename, bool defaultExp=false, const LibUtilities::FieldMetaDataMap &metadata=LibUtilities::NullFieldMetaDataMap)
void	Empty (int dim, int space)
void	FillGraph ()
void	FillBoundingBoxTree ()
std::vector< int >	GetElementsContainingPoint (PointGeomSharedPtr p)
void	ReadExpansionInfo ()
int	GetMeshDimension ()
	Dimension of the mesh (can be a 1D curve in 3D space). More...
int	GetSpaceDimension ()
	Dimension of the space (can be a 1D curve in 3D space). More...
void	SetDomainRange (NekDouble xmin, NekDouble xmax, NekDouble ymin=NekConstants::kNekUnsetDouble, NekDouble ymax=NekConstants::kNekUnsetDouble, NekDouble zmin=NekConstants::kNekUnsetDouble, NekDouble zmax=NekConstants::kNekUnsetDouble)
bool	CheckRange (Geometry2D &geom)
	Check if goemetry is in range definition if activated. More...
bool	CheckRange (Geometry3D &geom)
	Check if goemetry is in range definition if activated. More...
CompositeSharedPtr	GetComposite (int whichComposite)
GeometrySharedPtr	GetCompositeItem (int whichComposite, int whichItem)
void	GetCompositeList (const std::string &compositeStr, CompositeMap &compositeVector) const
std::map< int, CompositeSharedPtr > &	GetComposites ()
std::map< int, std::string > &	GetCompositesLabels ()
std::map< int, std::map< int, CompositeSharedPtr > > &	GetDomain ()
std::map< int, CompositeSharedPtr > &	GetDomain (int domain)
const ExpansionInfoMap &	GetExpansionInfo (const std::string variable="DefaultVar")
ExpansionInfoShPtr	GetExpansionInfo (GeometrySharedPtr geom, const std::string variable="DefaultVar")
void	SetExpansionInfo (std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef)
	Sets expansions given field definitions. More...
void	SetExpansionInfo (std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef, std::vector< std::vector< LibUtilities::PointsType >> &pointstype)
	Sets expansions given field definition, quadrature points. More...
void	SetExpansionInfoToEvenlySpacedPoints (int npoints=0)
	Sets expansions to have equispaced points. More...
void	SetExpansionInfoToNumModes (int nmodes)
	Reset expansion to have specified polynomial order nmodes. More...
void	SetExpansionInfoToPointOrder (int npts)
	Reset expansion to have specified point order npts. More...
void	SetExpansionInfo (const std::string variable, ExpansionInfoMapShPtr &exp)
	This function sets the expansion #exp in map with entry #variable. More...
void	SetSession (LibUtilities::SessionReaderSharedPtr pSession)
void	SetBasisKey (LibUtilities::ShapeType shape, LibUtilities::BasisKeyVector &keys, std::string var="DefaultVar")
	Sets the basis key for all expansions of the given shape. More...
void	ResetExpansionInfoToBasisKey (ExpansionInfoMapShPtr &expansionMap, LibUtilities::ShapeType shape, LibUtilities::BasisKeyVector &keys)
bool	SameExpansionInfo (const std::string var1, const std::string var2)
bool	ExpansionInfoDefined (const std::string var)
bool	CheckForGeomInfo (std::string parameter)
const std::string	GetGeomInfo (std::string parameter)
LibUtilities::BasisKeyVector	DefineBasisKeyFromExpansionTypeHomo (GeometrySharedPtr in, ExpansionType type_x, ExpansionType type_y, ExpansionType type_z, const int nummodes_x, const int nummodes_y, const int nummodes_z)
int	GetNvertices ()
PointGeomSharedPtr	GetVertex (int id)
SegGeomSharedPtr	GetSegGeom (int id)
CurveMap &	GetCurvedEdges ()
CurveMap &	GetCurvedFaces ()
std::map< int, PointGeomSharedPtr > &	GetAllPointGeoms ()
std::map< int, SegGeomSharedPtr > &	GetAllSegGeoms ()
TriGeomMap &	GetAllTriGeoms ()
QuadGeomMap &	GetAllQuadGeoms ()
TetGeomMap &	GetAllTetGeoms ()
PyrGeomMap &	GetAllPyrGeoms ()
PrismGeomMap &	GetAllPrismGeoms ()
HexGeomMap &	GetAllHexGeoms ()
int	GetNumElements ()
Geometry2DSharedPtr	GetGeometry2D (int gID)
LibUtilities::BasisKey	GetEdgeBasisKey (SegGeomSharedPtr edge, const std::string variable="DefaultVar")
GeometryLinkSharedPtr	GetElementsFromEdge (Geometry1DSharedPtr edge)
GeometryLinkSharedPtr	GetElementsFromFace (Geometry2DSharedPtr face)
LibUtilities::BasisKey	GetFaceBasisKey (Geometry2DSharedPtr face, const int facedir, const std::string variable="DefaultVar")
	3D functions More...
CompositeOrdering &	GetCompositeOrdering ()
void	SetCompositeOrdering (CompositeOrdering p_compOrder)
BndRegionOrdering &	GetBndRegionOrdering ()
void	SetBndRegionOrdering (BndRegionOrdering p_bndRegOrder)
void	ReadGeometry (LibUtilities::DomainRangeShPtr rng, bool fillGraph)
void	PartitionMesh (LibUtilities::SessionReaderSharedPtr session)
std::map< int, MeshEntity >	CreateMeshEntities ()
	Create mesh entities for this graph. More...
CompositeDescriptor	CreateCompositeDescriptor ()
MovementSharedPtr &	GetMovement ()

Static Public Member Functions
static MeshGraphSharedPtr	create ()
Static Public Member Functions inherited from Nektar::SpatialDomains::MeshGraph
static MeshGraphSharedPtr	Read (const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true)
static LibUtilities::BasisKeyVector	DefineBasisKeyFromExpansionType (GeometrySharedPtr in, ExpansionType type, const int order)

Static Public Attributes
static std::string	className

Protected Member Functions
virtual void	v_WriteGeometry (std::string &outfilename, bool defaultExp=false, const LibUtilities::FieldMetaDataMap &metadata=LibUtilities::NullFieldMetaDataMap) override
	Write out an XML file containing the GEOMETRY block representing this MeshGraph instance inside a NEKTAR tag. More...
virtual void	v_ReadGeometry (LibUtilities::DomainRangeShPtr rng, bool fillGraph) override
virtual void	v_PartitionMesh (LibUtilities::SessionReaderSharedPtr session) override
virtual void	v_ReadVertices ()
virtual void	v_ReadCurves ()
void	ReadDomain ()
virtual void	v_ReadEdges ()
virtual void	v_ReadFaces ()
void	ReadElements ()
void	ReadComposites ()
virtual void	v_ReadElements1D ()
virtual void	v_ReadElements2D ()
virtual void	v_ReadElements3D ()
void	ResolveGeomRef (const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
void	ResolveGeomRef1D (const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
void	ResolveGeomRef2D (const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
void	ResolveGeomRef3D (const std::string &prevToken, const std::string &token, CompositeSharedPtr &composite)
virtual void	v_WriteVertices (TiXmlElement *geomTag, PointGeomMap &verts)
virtual void	v_WriteEdges (TiXmlElement *geomTag, SegGeomMap &edges)
virtual void	v_WriteTris (TiXmlElement *faceTag, TriGeomMap &tris)
virtual void	v_WriteQuads (TiXmlElement *faceTag, QuadGeomMap &quads)
virtual void	v_WriteHexs (TiXmlElement *elmtTag, HexGeomMap &hexs)
virtual void	v_WritePrisms (TiXmlElement *elmtTag, PrismGeomMap &pris)
virtual void	v_WritePyrs (TiXmlElement *elmtTag, PyrGeomMap &pyrs)
virtual void	v_WriteTets (TiXmlElement *elmtTag, TetGeomMap &tets)
virtual void	v_WriteCurves (TiXmlElement *geomTag, CurveMap &edges, CurveMap &faces)
void	WriteComposites (TiXmlElement *geomTag, CompositeMap &comps, std::map< int, std::string > &compLabels)
void	WriteDomain (TiXmlElement *geomTag, std::map< int, CompositeMap > &domain)
void	WriteDefaultExpansion (TiXmlElement *root)
CompositeOrdering	CreateCompositeOrdering ()
Protected Member Functions inherited from Nektar::SpatialDomains::MeshGraph
void	PopulateFaceToElMap (Geometry3DSharedPtr element, int kNfaces)
	Given a 3D geometry object #element, populate the face to element map m_faceToElMap which maps faces to their corresponding element(s). More...
ExpansionInfoMapShPtr	SetUpExpansionInfoMap ()
std::string	GetCompositeString (CompositeSharedPtr comp)
	Returns a string representation of a composite. More...

Additional Inherited Members
Protected Attributes inherited from Nektar::SpatialDomains::MeshGraph
LibUtilities::SessionReaderSharedPtr	m_session
PointGeomMap	m_vertSet
CurveMap	m_curvedEdges
CurveMap	m_curvedFaces
SegGeomMap	m_segGeoms
TriGeomMap	m_triGeoms
QuadGeomMap	m_quadGeoms
TetGeomMap	m_tetGeoms
PyrGeomMap	m_pyrGeoms
PrismGeomMap	m_prismGeoms
HexGeomMap	m_hexGeoms
int	m_meshDimension
int	m_spaceDimension
int	m_partition
bool	m_meshPartitioned
CompositeMap	m_meshComposites
std::map< int, std::string >	m_compositesLabels
std::map< int, CompositeMap >	m_domain
LibUtilities::DomainRangeShPtr	m_domainRange
ExpansionInfoMapShPtrMap	m_expansionMapShPtrMap
GeomInfoMap	m_geomInfo
std::unordered_map< int, GeometryLinkSharedPtr >	m_faceToElMap
TiXmlElement *	m_xmlGeom
CompositeOrdering	m_compOrder
BndRegionOrdering	m_bndRegOrder
std::unique_ptr< GeomRTree >	m_boundingBoxTree
MovementSharedPtr	m_movement = nullptr

Detailed Description

Definition at line 47 of file MeshGraphXml.h.

Constructor & Destructor Documentation

MeshGraphXml()

Nektar::SpatialDomains::MeshGraphXml::MeshGraphXml ( )

inline

Definition at line 50 of file MeshGraphXml.h.

    {
    }

~MeshGraphXml()

virtual Nektar::SpatialDomains::MeshGraphXml::~MeshGraphXml ( )

inlinevirtual

Definition at line 54 of file MeshGraphXml.h.

    {
    }

Member Function Documentation

create()

static MeshGraphSharedPtr Nektar::SpatialDomains::MeshGraphXml::create ( )

inlinestatic

Definition at line 62 of file MeshGraphXml.h.

    {
        return MemoryManager<MeshGraphXml>::AllocateSharedPtr();
    }

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

CreateCompositeOrdering()

CompositeOrdering Nektar::SpatialDomains::MeshGraphXml::CreateCompositeOrdering ( )

protected

Definition at line 3231 of file MeshGraphXml.cpp.

{
    CompositeOrdering ret;
 
    for (auto &c : m_meshComposites)
    {
        bool fillComp = true;
        for (auto &d : m_domain[0])
        {
            if (c.second == d.second)
            {
                fillComp = false;
            }
        }
        if (fillComp)
        {
            std::vector<unsigned int> ids;
            for (auto &elmt : c.second->m_geomVec)
            {
                ids.push_back(elmt->GetGlobalID());
            }
            ret[c.first] = ids;
        }
    }
 
    return ret;
}

ReadComposites()

void Nektar::SpatialDomains::MeshGraphXml::ReadComposites ( )

protected

Look for elements in ELEMENT block.

All elements are of the form: "<? ID="#"> ... </?>", with ? being the element type.

Parse out the element components corresponding to type of element.

Keep looking for additional composite definitions.

Definition at line 1854 of file MeshGraphXml.cpp.

{
    TiXmlElement *field = NULL;
 
    /// Look for elements in ELEMENT block.
    field = m_xmlGeom->FirstChildElement("COMPOSITE");
 
    ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
 
    TiXmlElement *node = field->FirstChildElement("C");
 
    // Sequential counter for the composite numbers.
    int nextCompositeNumber = -1;
 
    while (node)
    {
        /// All elements are of the form: "<? ID="#"> ... </?>", with
        /// ? being the element type.
 
        nextCompositeNumber++;
 
        int indx;
        int err = node->QueryIntAttribute("ID", &indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
        // ASSERTL0(indx == nextCompositeNumber, "Composite IDs must begin with
        // zero and be sequential.");
 
        TiXmlNode *compositeChild = node->FirstChild();
        // This is primarily to skip comments that may be present.
        // Comments appear as nodes just like elements.
        // We are specifically looking for text in the body
        // of the definition.
        while (compositeChild &&
               compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
        {
            compositeChild = compositeChild->NextSibling();
        }
 
        ASSERTL0(compositeChild, "Unable to read composite definition body.");
        std::string compositeStr = compositeChild->ToText()->ValueStr();
 
        /// Parse out the element components corresponding to type of element.
 
        std::istringstream compositeDataStrm(compositeStr.c_str());
 
        try
        {
            bool first = true;
            std::string prevCompositeElementStr;
 
            while (!compositeDataStrm.fail())
            {
                std::string compositeElementStr;
                compositeDataStrm >> compositeElementStr;
 
                if (!compositeDataStrm.fail())
                {
                    if (first)
                    {
                        first = false;
 
                        CompositeSharedPtr curVector =
                            MemoryManager<Composite>::AllocateSharedPtr();
                        m_meshComposites[indx] = curVector;
                    }
 
                    if (compositeElementStr.length() > 0)
                    {
                        ResolveGeomRef(prevCompositeElementStr,
                                       compositeElementStr,
                                       m_meshComposites[indx]);
                    }
                    prevCompositeElementStr = compositeElementStr;
                }
            }
        }
        catch (...)
        {
            NEKERROR(
                ErrorUtil::efatal,
                (std::string("Unable to read COMPOSITE data for composite: ") +
                 compositeStr)
                    .c_str());
        }
 
        // Read optional name as string and save to m_compositeLabels if exists
        std::string name;
        err = node->QueryStringAttribute("NAME", &name);
        if (err == TIXML_SUCCESS)
        {
            m_compositesLabels[indx] = name;
        }
 
        /// Keep looking for additional composite definitions.
        node = node->NextSiblingElement("C");
    }
 
    ASSERTL0(nextCompositeNumber >= 0,
             "At least one composite must be specified.");
}

References ASSERTL0, CellMLToNektar.pycml::name, and NEKERROR.

ReadDomain()

void Nektar::SpatialDomains::MeshGraphXml::ReadDomain ( )

protected

Look for data in DOMAIN block.

Elements are of the form: "<D ID = "N"> ... </D>". Read the ID field first.

Keep looking

Definition at line 945 of file MeshGraphXml.cpp.

{
    TiXmlElement *domain = NULL;
    /// Look for data in DOMAIN block.
    domain = m_xmlGeom->FirstChildElement("DOMAIN");
 
    ASSERTL0(domain, "Unable to find DOMAIN tag in file.");
 
    /// Elements are of the form: "<D ID = "N"> ... </D>".
    /// Read the ID field first.
    TiXmlElement *multidomains = domain->FirstChildElement("D");
 
    if (multidomains)
    {
        while (multidomains)
        {
            int indx;
            int err = multidomains->QueryIntAttribute("ID", &indx);
            ASSERTL0(err == TIXML_SUCCESS,
                     "Unable to read attribute ID in Domain.");
 
            TiXmlNode *elementChild = multidomains->FirstChild();
            while (elementChild &&
                   elementChild->Type() != TiXmlNode::TINYXML_TEXT)
            {
                elementChild = elementChild->NextSibling();
            }
 
            ASSERTL0(elementChild, "Unable to read DOMAIN body.");
            std::string elementStr = elementChild->ToText()->ValueStr();
 
            elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
 
            std::string::size_type indxBeg = elementStr.find_first_of('[') + 1;
            std::string::size_type indxEnd = elementStr.find_last_of(']') - 1;
            std::string indxStr =
                elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
 
            ASSERTL0(
                !indxStr.empty(),
                "Unable to read domain's composite index (index missing?).");
 
            // Read the domain composites.
            // Parse the composites into a list.
            map<int, CompositeSharedPtr> unrollDomain;
            GetCompositeList(indxStr, unrollDomain);
            m_domain[indx] = unrollDomain;
 
            ASSERTL0(!m_domain[indx].empty(),
                     (std::string(
                          "Unable to obtain domain's referenced composite: ") +
                      indxStr)
                         .c_str());
 
            /// Keep looking
            multidomains = multidomains->NextSiblingElement("D");
        }
    }
    else // previous definition of just one composite
    {
 
        // find the non comment portion of the body.
        TiXmlNode *elementChild = domain->FirstChild();
        while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
        {
            elementChild = elementChild->NextSibling();
        }
 
        ASSERTL0(elementChild, "Unable to read DOMAIN body.");
        std::string elementStr = elementChild->ToText()->ValueStr();
 
        elementStr = elementStr.substr(elementStr.find_first_not_of(" "));
 
        std::string::size_type indxBeg = elementStr.find_first_of('[') + 1;
        std::string::size_type indxEnd = elementStr.find_last_of(']') - 1;
        std::string indxStr = elementStr.substr(indxBeg, indxEnd - indxBeg + 1);
 
        ASSERTL0(!indxStr.empty(),
                 "Unable to read domain's composite index (index missing?).");
 
        // Read the domain composites.
        // Parse the composites into a list.
        map<int, CompositeSharedPtr> fullDomain;
        GetCompositeList(indxStr, fullDomain);
        m_domain[0] = fullDomain;
 
        ASSERTL0(
            !m_domain[0].empty(),
            (std::string("Unable to obtain domain's referenced composite: ") +
             indxStr)
                .c_str());
    }
}

References ASSERTL0.

ReadElements()

void Nektar::SpatialDomains::MeshGraphXml::ReadElements ( )

protected

Definition at line 1267 of file MeshGraphXml.cpp.

{
    switch (m_meshDimension)
    {
        case 1:
            v_ReadElements1D();
            break;
        case 2:
            v_ReadElements2D();
            break;
        case 3:
            v_ReadElements3D();
            break;
    }
}

ResolveGeomRef()

void Nektar::SpatialDomains::MeshGraphXml::ResolveGeomRef ( const std::string &  prevToken, const std::string &  token, CompositeSharedPtr &  composite  )

protected

Definition at line 1955 of file MeshGraphXml.cpp.

{
    switch (m_meshDimension)
    {
        case 1:
            ResolveGeomRef1D(prevToken, token, composite);
            break;
        case 2:
            ResolveGeomRef2D(prevToken, token, composite);
            break;
        case 3:
            ResolveGeomRef3D(prevToken, token, composite);
            break;
    }
}

ResolveGeomRef1D()

void Nektar::SpatialDomains::MeshGraphXml::ResolveGeomRef1D ( const std::string &  prevToken, const std::string &  token, CompositeSharedPtr &  composite  )

protected

Definition at line 1973 of file MeshGraphXml.cpp.

{
    try
    {
        std::istringstream tokenStream(token);
        std::istringstream prevTokenStream(prevToken);
 
        char type;
        char prevType;
 
        tokenStream >> type;
 
        std::string::size_type indxBeg = token.find_first_of('[') + 1;
        std::string::size_type indxEnd = token.find_last_of(']') - 1;
 
        ASSERTL0(
            indxBeg <= indxEnd,
            (std::string("Error reading index definition:") + token).c_str());
 
        std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
 
        typedef vector<unsigned int> SeqVectorType;
        SeqVectorType seqVector;
 
        if (!ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector))
        {
            NEKERROR(ErrorUtil::efatal,
                     (std::string("Ill-formed sequence definition: ") + indxStr)
                         .c_str());
        }
 
        prevTokenStream >> prevType;
 
        // All composites must be of the same dimension.
        bool validSequence =
            (prevToken.empty() || // No previous, then current is just fine.
             (type == 'V' && prevType == 'V') ||
             (type == 'S' && prevType == 'S'));
 
        ASSERTL0(validSequence,
                 std::string("Invalid combination of composite items: ") +
                     type + " and " + prevType + ".");
 
        switch (type)
        {
            case 'V': // Vertex
                for (SeqVectorType::iterator iter = seqVector.begin();
                     iter != seqVector.end(); ++iter)
                {
                    if (m_vertSet.find(*iter) == m_vertSet.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown vertex index: " +
                                     std::to_string(*iter));
                    }
                    else
                    {
                        composite->m_geomVec.push_back(m_vertSet[*iter]);
                    }
                }
                break;
 
            case 'S': // Segment
                for (SeqVectorType::iterator iter = seqVector.begin();
                     iter != seqVector.end(); ++iter)
                {
                    if (m_segGeoms.find(*iter) == m_segGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown segment index: " +
                                     std::to_string(*iter));
                    }
                    else
                    {
                        composite->m_geomVec.push_back(m_segGeoms[*iter]);
                    }
                }
                break;
 
            default:
                NEKERROR(ErrorUtil::efatal,
                         "Unrecognized composite token: " + token);
        }
    }
    catch (...)
    {
        NEKERROR(ErrorUtil::efatal,
                 "Problem processing composite token: " + token);
    }
 
    return;
}

References ASSERTL0, and NEKERROR.

ResolveGeomRef2D()

void Nektar::SpatialDomains::MeshGraphXml::ResolveGeomRef2D ( const std::string &  prevToken, const std::string &  token, CompositeSharedPtr &  composite  )

protected

Definition at line 2068 of file MeshGraphXml.cpp.

{
    try
    {
        std::istringstream tokenStream(token);
        std::istringstream prevTokenStream(prevToken);
 
        char type;
        char prevType;
 
        tokenStream >> type;
 
        std::string::size_type indxBeg = token.find_first_of('[') + 1;
        std::string::size_type indxEnd = token.find_last_of(']') - 1;
 
        ASSERTL0(
            indxBeg <= indxEnd,
            (std::string("Error reading index definition:") + token).c_str());
 
        std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
        std::vector<unsigned int> seqVector;
        std::vector<unsigned int>::iterator seqIter;
 
        bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
 
        ASSERTL0(err,
                 (std::string("Error reading composite elements: ") + indxStr)
                     .c_str());
 
        prevTokenStream >> prevType;
 
        // All composites must be of the same dimension.
        bool validSequence =
            (prevToken.empty() || // No previous, then current is just fine.
             (type == 'V' && prevType == 'V') ||
             (type == 'E' && prevType == 'E') ||
             ((type == 'T' || type == 'Q') &&
              (prevType == 'T' || prevType == 'Q')));
 
        ASSERTL0(validSequence,
                 std::string("Invalid combination of composite items: ") +
                     type + " and " + prevType + ".");
 
        switch (type)
        {
            case 'E': // Edge
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown edge index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        composite->m_geomVec.push_back(m_segGeoms[*seqIter]);
                    }
                }
                break;
 
            case 'T': // Triangle
            {
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_triGeoms.count(*seqIter) == 0)
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown triangle index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*m_triGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_triGeoms[*seqIter]);
                        }
                    }
                }
            }
            break;
 
            case 'Q': // Quad
            {
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_quadGeoms.count(*seqIter) == 0)
                    {
                        NEKERROR(
                            ErrorUtil::ewarning,
                            "Unknown quad index: " + std::to_string(*seqIter) +
                                " in Composite section");
                    }
                    else
                    {
                        if (CheckRange(*m_quadGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_quadGeoms[*seqIter]);
                        }
                    }
                }
            }
            break;
 
            case 'V': // Vertex
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (*seqIter >= m_vertSet.size())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown vertex index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        composite->m_geomVec.push_back(m_vertSet[*seqIter]);
                    }
                }
                break;
 
            default:
                NEKERROR(ErrorUtil::efatal,
                         "Unrecognized composite token: " + token);
        }
    }
    catch (...)
    {
        NEKERROR(ErrorUtil::efatal,
                 "Problem processing composite token: " + token);
    }
 
    return;
}

References ASSERTL0, and NEKERROR.

ResolveGeomRef3D()

void Nektar::SpatialDomains::MeshGraphXml::ResolveGeomRef3D ( const std::string &  prevToken, const std::string &  token, CompositeSharedPtr &  composite  )

protected

Definition at line 2210 of file MeshGraphXml.cpp.

{
    try
    {
        std::istringstream tokenStream(token);
        std::istringstream prevTokenStream(prevToken);
 
        char type;
        char prevType;
 
        tokenStream >> type;
 
        std::string::size_type indxBeg = token.find_first_of('[') + 1;
        std::string::size_type indxEnd = token.find_last_of(']') - 1;
 
        ASSERTL0(indxBeg <= indxEnd,
                 "Error reading index definition: " + token);
 
        std::string indxStr = token.substr(indxBeg, indxEnd - indxBeg + 1);
 
        std::vector<unsigned int> seqVector;
        std::vector<unsigned int>::iterator seqIter;
 
        bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
 
        ASSERTL0(err, "Error reading composite elements: " + indxStr);
 
        prevTokenStream >> prevType;
 
        // All composites must be of the same dimension.  This map makes things
        // clean to compare.
        map<char, int> typeMap;
        typeMap['V'] = 1; // Vertex
        typeMap['E'] = 1; // Edge
        typeMap['T'] = 2; // Triangle
        typeMap['Q'] = 2; // Quad
        typeMap['A'] = 3; // Tet
        typeMap['P'] = 3; // Pyramid
        typeMap['R'] = 3; // Prism
        typeMap['H'] = 3; // Hex
 
        // Make sure only geoms of the same dimension are combined.
        bool validSequence =
            (prevToken.empty() || (typeMap[type] == typeMap[prevType]));
 
        ASSERTL0(validSequence,
                 std::string("Invalid combination of composite items: ") +
                     type + " and " + prevType + ".");
 
        switch (type)
        {
            case 'V': // Vertex
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_vertSet.find(*seqIter) == m_vertSet.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown vertex index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        composite->m_geomVec.push_back(m_vertSet[*seqIter]);
                    }
                }
                break;
 
            case 'E': // Edge
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_segGeoms.find(*seqIter) == m_segGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown edge index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        composite->m_geomVec.push_back(m_segGeoms[*seqIter]);
                    }
                }
                break;
 
            case 'F': // Face
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    Geometry2DSharedPtr face = GetGeometry2D(*seqIter);
                    if (face == Geometry2DSharedPtr())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown face index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*face))
                        {
                            composite->m_geomVec.push_back(face);
                        }
                    }
                }
                break;
 
            case 'T': // Triangle
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_triGeoms.find(*seqIter) == m_triGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown triangle index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*m_triGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_triGeoms[*seqIter]);
                        }
                    }
                }
                break;
 
            case 'Q': // Quad
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_quadGeoms.find(*seqIter) == m_quadGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown quad index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*m_quadGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_quadGeoms[*seqIter]);
                        }
                    }
                }
                break;
 
            // Tetrahedron
            case 'A':
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_tetGeoms.find(*seqIter) == m_tetGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown tet index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*m_tetGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_tetGeoms[*seqIter]);
                        }
                    }
                }
                break;
 
            // Pyramid
            case 'P':
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_pyrGeoms.find(*seqIter) == m_pyrGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown pyramid index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*m_pyrGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_pyrGeoms[*seqIter]);
                        }
                    }
                }
                break;
 
            // Prism
            case 'R':
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_prismGeoms.find(*seqIter) == m_prismGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown prism index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*m_prismGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_prismGeoms[*seqIter]);
                        }
                    }
                }
                break;
 
            // Hex
            case 'H':
                for (seqIter = seqVector.begin(); seqIter != seqVector.end();
                     ++seqIter)
                {
                    if (m_hexGeoms.find(*seqIter) == m_hexGeoms.end())
                    {
                        NEKERROR(ErrorUtil::ewarning,
                                 "Unknown hex index: " +
                                     std::to_string(*seqIter));
                    }
                    else
                    {
                        if (CheckRange(*m_hexGeoms[*seqIter]))
                        {
                            composite->m_geomVec.push_back(
                                m_hexGeoms[*seqIter]);
                        }
                    }
                }
                break;
 
            default:
                NEKERROR(ErrorUtil::efatal,
                         "Unrecognized composite token: " + token);
        }
    }
    catch (...)
    {
        NEKERROR(ErrorUtil::efatal,
                 "Problem processing composite token: " + token);
    }
 
    return;
}

References ASSERTL0, and NEKERROR.

v_PartitionMesh()

void Nektar::SpatialDomains::MeshGraphXml::v_PartitionMesh ( LibUtilities::SessionReaderSharedPtr  session )

overrideprotectedvirtual

Implements Nektar::SpatialDomains::MeshGraph.

Definition at line 61 of file MeshGraphXml.cpp.

{
    // Get row of comm, or the whole comm if not split
    LibUtilities::CommSharedPtr comm     = session->GetComm();
    LibUtilities::CommSharedPtr commMesh = comm->GetRowComm();
    const bool isRoot                    = comm->TreatAsRankZero();
 
    m_session = session;
 
    // Load file for root process only (since this is always needed)
    // and determine if the provided geometry has already been
    // partitioned. This will be the case if the user provides the
    // directory of mesh partitions as an input. Partitioned geometries
    // have the attribute
    //    PARTITION=X
    // where X is the number of the partition (and should match the
    // process rank). The result is shared with all other processes.
    int isPartitioned = 0;
    if (isRoot)
    {
        if (m_session->DefinesElement("Nektar/Geometry"))
        {
            if (m_session->GetElement("Nektar/Geometry")
                    ->Attribute("PARTITION"))
            {
                std::cout << "Using pre-partitioned mesh." << std::endl;
                isPartitioned = 1;
            }
        }
    }
    comm->Bcast(isPartitioned, 0);
 
    // If the mesh is already partitioned, we are done. Remaining
    // processes must load their partitions.
    if (isPartitioned)
    {
        if (!isRoot)
        {
            m_session->InitSession();
        }
    }
    else
    {
        // Default partitioner to use is Metis. Use Scotch as default if it is
        // installed. Override default with command-line flags if they are set.
        string partitionerName = "Metis";
        if (GetMeshPartitionFactory().ModuleExists("Scotch"))
        {
            partitionerName = "Scotch";
        }
        if (session->DefinesCmdLineArgument("use-metis"))
        {
            partitionerName = "Metis";
        }
        if (session->DefinesCmdLineArgument("use-scotch"))
        {
            partitionerName = "Scotch";
        }
 
        // Mesh has not been partitioned so do partitioning if required.  Note
        // in the serial case nothing is done as we have already loaded the
        // mesh.
        if (session->DefinesCmdLineArgument("part-only") ||
            session->DefinesCmdLineArgument("part-only-overlapping"))
        {
            // Perform partitioning of the mesh only. For this we insist the
            // code is run in serial (parallel execution is pointless).
            ASSERTL0(comm->GetSize() == 1,
                     "The 'part-only' option should be used in serial.");
 
            // Read 'lite' geometry information
            ReadGeometry(LibUtilities::NullDomainRangeShPtr, false);
 
            // Number of partitions is specified by the parameter.
            int nParts;
            auto comp = CreateCompositeDescriptor();
 
            MeshPartitionSharedPtr partitioner =
                GetMeshPartitionFactory().CreateInstance(
                    partitionerName, session, session->GetComm(),
                    m_meshDimension, CreateMeshEntities(), comp);
 
            if (session->DefinesCmdLineArgument("part-only"))
            {
                nParts = session->GetCmdLineArgument<int>("part-only");
                partitioner->PartitionMesh(nParts, true);
            }
            else
            {
                nParts =
                    session->GetCmdLineArgument<int>("part-only-overlapping");
                partitioner->PartitionMesh(nParts, true, true);
            }
 
            vector<set<unsigned int>> elmtIDs;
            vector<unsigned int> parts(nParts);
            for (int i = 0; i < nParts; ++i)
            {
                vector<unsigned int> elIDs;
                set<unsigned int> tmp;
                partitioner->GetElementIDs(i, elIDs);
                tmp.insert(elIDs.begin(), elIDs.end());
                elmtIDs.push_back(tmp);
                parts[i] = i;
            }
 
            this->WriteXMLGeometry(m_session->GetSessionName(), elmtIDs, parts);
 
            if (isRoot && session->DefinesCmdLineArgument("part-info"))
            {
                partitioner->PrintPartInfo(std::cout);
            }
 
            session->Finalise();
            exit(0);
        }
 
        if (commMesh->GetSize() > 1)
        {
            int nParts = commMesh->GetSize();
 
            if (session->GetSharedFilesystem())
            {
                vector<unsigned int> keys, vals;
                int i;
 
                if (isRoot)
                {
                    // Read 'lite' geometry information
                    ReadGeometry(LibUtilities::NullDomainRangeShPtr, false);
 
                    // Store composite ordering and boundary information.
                    m_compOrder = CreateCompositeOrdering();
                    auto comp   = CreateCompositeDescriptor();
 
                    // Create mesh partitioner.
                    MeshPartitionSharedPtr partitioner =
                        GetMeshPartitionFactory().CreateInstance(
                            partitionerName, session, session->GetComm(),
                            m_meshDimension, CreateMeshEntities(), comp);
 
                    partitioner->PartitionMesh(nParts, true);
 
                    vector<set<unsigned int>> elmtIDs;
                    vector<unsigned int> parts(nParts);
                    for (i = 0; i < nParts; ++i)
                    {
                        vector<unsigned int> elIDs;
                        set<unsigned int> tmp;
                        partitioner->GetElementIDs(i, elIDs);
                        tmp.insert(elIDs.begin(), elIDs.end());
                        elmtIDs.push_back(tmp);
                        parts[i] = i;
                    }
 
                    // Call WriteGeometry to write out partition files. This
                    // will populate m_bndRegOrder.
                    this->WriteXMLGeometry(m_session->GetSessionName(), elmtIDs,
                                           parts);
 
                    // Communicate orderings to the other processors.
 
                    // First send sizes of the orderings and boundary
                    // regions to allocate storage on the remote end.
                    keys.resize(2);
                    keys[0] = m_compOrder.size();
                    keys[1] = m_bndRegOrder.size();
                    comm->Bcast(keys, 0);
 
                    // Construct the keys and sizes of values for composite
                    // ordering
                    keys.resize(m_compOrder.size());
                    vals.resize(m_compOrder.size());
 
                    i = 0;
                    for (auto &cIt : m_compOrder)
                    {
                        keys[i]   = cIt.first;
                        vals[i++] = cIt.second.size();
                    }
 
                    // Send across data.
                    comm->Bcast(keys, 0);
                    comm->Bcast(vals, 0);
                    for (auto &cIt : m_compOrder)
                    {
                        comm->Bcast(cIt.second, 0);
                    }
 
                    // Construct the keys and sizes of values for composite
                    // ordering
                    keys.resize(m_bndRegOrder.size());
                    vals.resize(m_bndRegOrder.size());
 
                    i = 0;
                    for (auto &bIt : m_bndRegOrder)
                    {
                        keys[i]   = bIt.first;
                        vals[i++] = bIt.second.size();
                    }
 
                    // Send across data.
                    if (!keys.empty())
                    {
                        comm->Bcast(keys, 0);
                    }
                    if (!vals.empty())
                    {
                        comm->Bcast(vals, 0);
                    }
                    for (auto &bIt : m_bndRegOrder)
                    {
                        comm->Bcast(bIt.second, 0);
                    }
 
                    if (session->DefinesCmdLineArgument("part-info"))
                    {
                        partitioner->PrintPartInfo(std::cout);
                    }
                }
                else
                {
                    keys.resize(2);
                    comm->Bcast(keys, 0);
 
                    int cmpSize = keys[0];
                    int bndSize = keys[1];
 
                    keys.resize(cmpSize);
                    vals.resize(cmpSize);
                    comm->Bcast(keys, 0);
                    comm->Bcast(vals, 0);
 
                    for (int i = 0; i < keys.size(); ++i)
                    {
                        vector<unsigned int> tmp(vals[i]);
                        comm->Bcast(tmp, 0);
                        m_compOrder[keys[i]] = tmp;
                    }
 
                    keys.resize(bndSize);
                    vals.resize(bndSize);
                    if (!keys.empty())
                    {
                        comm->Bcast(keys, 0);
                    }
                    if (!vals.empty())
                    {
                        comm->Bcast(vals, 0);
                    }
                    for (int i = 0; i < keys.size(); ++i)
                    {
                        vector<unsigned int> tmp(vals[i]);
                        comm->Bcast(tmp, 0);
                        m_bndRegOrder[keys[i]] = tmp;
                    }
                }
            }
            else
            {
                m_session->InitSession();
                ReadGeometry(LibUtilities::NullDomainRangeShPtr, false);
 
                m_compOrder = CreateCompositeOrdering();
                auto comp   = CreateCompositeDescriptor();
 
                // Partitioner now operates in parallel. Each process receives
                // partitioning over interconnect and writes its own session
                // file to the working directory.
                MeshPartitionSharedPtr partitioner =
                    GetMeshPartitionFactory().CreateInstance(
                        partitionerName, session, session->GetComm(),
                        m_meshDimension, CreateMeshEntities(), comp);
 
                partitioner->PartitionMesh(nParts, false);
 
                vector<unsigned int> parts(1), tmp;
                parts[0] = commMesh->GetRank();
                vector<set<unsigned int>> elIDs(1);
                partitioner->GetElementIDs(parts[0], tmp);
                elIDs[0].insert(tmp.begin(), tmp.end());
 
                this->WriteXMLGeometry(session->GetSessionName(), elIDs, parts);
 
                if (m_session->DefinesCmdLineArgument("part-info") && isRoot)
                {
                    partitioner->PrintPartInfo(std::cout);
                }
            }
 
            // Wait for all processors to finish their writing activities.
            comm->Block();
 
            std::string dirname = m_session->GetSessionName() + "_xml";
            fs::path pdirname(dirname);
            boost::format pad("P%1$07d.xml");
            pad % comm->GetRowComm()->GetRank();
            fs::path pFilename(pad.str());
            fs::path fullpath = pdirname / pFilename;
 
            std::vector<std::string> filenames = {
                LibUtilities::PortablePath(fullpath)};
            m_session->InitSession(filenames);
        }
        else if (!isRoot)
        {
            // No partitioning, non-root processors need to read the session
            // file -- handles case where --npz is the same as number of
            // processors.
            m_session->InitSession();
        }
    }
}

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), CellMLToNektar.pycml::format, Nektar::SpatialDomains::GetMeshPartitionFactory(), Nektar::LibUtilities::NullDomainRangeShPtr, and Nektar::LibUtilities::PortablePath().

v_ReadCurves()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadCurves ( )

protectedvirtual

Look for elements in CURVE block.

All curves are of the form: "<? ID="#" TYPE="GLL OR other points type" NUMPOINTS="#"> ... </?>", with ? being an element type (either E or F).

These should be ordered.

Read id attribute.

Read edge id attribute.

Read text edgelement description.

Parse out the element components corresponding to type of element.

Read id attribute.

Read face id attribute.

Read text face element description.

Parse out the element components corresponding to type of element.

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 625 of file MeshGraphXml.cpp.

{
    // check to see if any scaling parameters are in
    // attributes and determine these values
    TiXmlElement *element = m_xmlGeom->FirstChildElement("VERTEX");
    ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
 
    NekDouble xscale, yscale, zscale;
 
    LibUtilities::Interpreter expEvaluator;
    const char *xscal = element->Attribute("XSCALE");
    if (!xscal)
    {
        xscale = 1.0;
    }
    else
    {
        std::string xscalstr = xscal;
        int expr_id          = expEvaluator.DefineFunction("", xscalstr);
        xscale               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *yscal = element->Attribute("YSCALE");
    if (!yscal)
    {
        yscale = 1.0;
    }
    else
    {
        std::string yscalstr = yscal;
        int expr_id          = expEvaluator.DefineFunction("", yscalstr);
        yscale               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *zscal = element->Attribute("ZSCALE");
    if (!zscal)
    {
        zscale = 1.0;
    }
    else
    {
        std::string zscalstr = zscal;
        int expr_id          = expEvaluator.DefineFunction("", zscalstr);
        zscale               = expEvaluator.Evaluate(expr_id);
    }
 
    NekDouble xmove, ymove, zmove;
 
    // check to see if any moving parameters are in
    // attributes and determine these values
 
    const char *xmov = element->Attribute("XMOVE");
    if (!xmov)
    {
        xmove = 0.0;
    }
    else
    {
        std::string xmovstr = xmov;
        int expr_id         = expEvaluator.DefineFunction("", xmovstr);
        xmove               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *ymov = element->Attribute("YMOVE");
    if (!ymov)
    {
        ymove = 0.0;
    }
    else
    {
        std::string ymovstr = ymov;
        int expr_id         = expEvaluator.DefineFunction("", ymovstr);
        ymove               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *zmov = element->Attribute("ZMOVE");
    if (!zmov)
    {
        zmove = 0.0;
    }
    else
    {
        std::string zmovstr = zmov;
        int expr_id         = expEvaluator.DefineFunction("", zmovstr);
        zmove               = expEvaluator.Evaluate(expr_id);
    }
 
    int err;
 
    /// Look for elements in CURVE block.
    TiXmlElement *field = m_xmlGeom->FirstChildElement("CURVED");
 
    if (!field) // return if no curved entities
    {
        return;
    }
 
    /// All curves are of the form: "<? ID="#" TYPE="GLL OR other
    /// points type" NUMPOINTS="#"> ... </?>", with ? being an
    /// element type (either E or F).
 
    TiXmlElement *edgelement = field->FirstChildElement("E");
 
    int edgeindx, edgeid;
    int nextEdgeNumber = -1;
 
    while (edgelement)
    {
        /// These should be ordered.
        nextEdgeNumber++;
 
        std::string edge(edgelement->ValueStr());
        ASSERTL0(edge == "E",
                 (std::string("Unknown 3D curve type:") + edge).c_str());
 
        /// Read id attribute.
        err = edgelement->QueryIntAttribute("ID", &edgeindx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute ID.");
 
        /// Read edge id attribute.
        err = edgelement->QueryIntAttribute("EDGEID", &edgeid);
        ASSERTL0(err == TIXML_SUCCESS,
                 "Unable to read curve attribute EDGEID.");
 
        /// Read text edgelement description.
        std::string elementStr;
        TiXmlNode *elementChild = edgelement->FirstChild();
 
        while (elementChild)
        {
            // Accumulate all non-comment element data
            if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
            {
                elementStr += elementChild->ToText()->ValueStr();
                elementStr += " ";
            }
            elementChild = elementChild->NextSibling();
        }
 
        ASSERTL0(!elementStr.empty(), "Unable to read curve description body.");
 
        /// Parse out the element components corresponding to type of
        /// element.
        if (edge == "E")
        {
            int numPts = 0;
            // Determine the points type
            std::string typeStr = edgelement->Attribute("TYPE");
            ASSERTL0(!typeStr.empty(), "TYPE must be specified in "
                                       "points definition");
 
            LibUtilities::PointsType type;
            const std::string *begStr = LibUtilities::kPointsTypeStr;
            const std::string *endStr =
                LibUtilities::kPointsTypeStr + LibUtilities::SIZE_PointsType;
            const std::string *ptsStr = std::find(begStr, endStr, typeStr);
 
            ASSERTL0(ptsStr != endStr, "Invalid points type.");
            type = (LibUtilities::PointsType)(ptsStr - begStr);
 
            // Determine the number of points
            err = edgelement->QueryIntAttribute("NUMPOINTS", &numPts);
            ASSERTL0(err == TIXML_SUCCESS,
                     "Unable to read curve attribute NUMPOINTS.");
            CurveSharedPtr curve(
                MemoryManager<Curve>::AllocateSharedPtr(edgeid, type));
 
            // Read points (x, y, z)
            NekDouble xval, yval, zval;
            std::istringstream elementDataStrm(elementStr.c_str());
            try
            {
                while (!elementDataStrm.fail())
                {
                    elementDataStrm >> xval >> yval >> zval;
 
                    xval = xval * xscale + xmove;
                    yval = yval * yscale + ymove;
                    zval = zval * zscale + zmove;
 
                    // Need to check it here because we may not be
                    // good after the read indicating that there
                    // was nothing to read.
                    if (!elementDataStrm.fail())
                    {
                        PointGeomSharedPtr vert(
                            MemoryManager<PointGeom>::AllocateSharedPtr(
                                m_meshDimension, edgeindx, xval, yval, zval));
 
                        curve->m_points.push_back(vert);
                    }
                }
            }
            catch (...)
            {
                NEKERROR(ErrorUtil::efatal,
                         (std::string("Unable to read curve data for EDGE: ") +
                          elementStr)
                             .c_str());
            }
 
            ASSERTL0(curve->m_points.size() == numPts,
                     "Number of points specificed by attribute "
                     "NUMPOINTS is different from number of points "
                     "in list (edgeid = " +
                         boost::lexical_cast<string>(edgeid));
 
            m_curvedEdges[edgeid] = curve;
 
            edgelement = edgelement->NextSiblingElement("E");
 
        } // end if-loop
 
    } // end while-loop
 
    TiXmlElement *facelement = field->FirstChildElement("F");
    int faceindx, faceid;
 
    while (facelement)
    {
        std::string face(facelement->ValueStr());
        ASSERTL0(face == "F",
                 (std::string("Unknown 3D curve type: ") + face).c_str());
 
        /// Read id attribute.
        err = facelement->QueryIntAttribute("ID", &faceindx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read curve attribute ID.");
 
        /// Read face id attribute.
        err = facelement->QueryIntAttribute("FACEID", &faceid);
        ASSERTL0(err == TIXML_SUCCESS,
                 "Unable to read curve attribute FACEID.");
 
        /// Read text face element description.
        std::string elementStr;
        TiXmlNode *elementChild = facelement->FirstChild();
 
        while (elementChild)
        {
            // Accumulate all non-comment element data
            if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
            {
                elementStr += elementChild->ToText()->ValueStr();
                elementStr += " ";
            }
            elementChild = elementChild->NextSibling();
        }
 
        ASSERTL0(!elementStr.empty(), "Unable to read curve description body.");
 
        /// Parse out the element components corresponding to type of
        /// element.
        if (face == "F")
        {
            std::string typeStr = facelement->Attribute("TYPE");
            ASSERTL0(!typeStr.empty(), "TYPE must be specified in "
                                       "points definition");
            LibUtilities::PointsType type;
            const std::string *begStr = LibUtilities::kPointsTypeStr;
            const std::string *endStr =
                LibUtilities::kPointsTypeStr + LibUtilities::SIZE_PointsType;
            const std::string *ptsStr = std::find(begStr, endStr, typeStr);
 
            ASSERTL0(ptsStr != endStr, "Invalid points type.");
            type = (LibUtilities::PointsType)(ptsStr - begStr);
 
            std::string numptsStr = facelement->Attribute("NUMPOINTS");
            ASSERTL0(!numptsStr.empty(),
                     "NUMPOINTS must be specified in points definition");
            int numPts = 0;
            std::stringstream s;
            s << numptsStr;
            s >> numPts;
 
            CurveSharedPtr curve(
                MemoryManager<Curve>::AllocateSharedPtr(faceid, type));
 
            ASSERTL0(numPts >= 3, "NUMPOINTS for face must be greater than 2");
 
            if (numPts == 3)
            {
                ASSERTL0(ptsStr != endStr, "Invalid points type.");
            }
 
            // Read points (x, y, z)
            NekDouble xval, yval, zval;
            std::istringstream elementDataStrm(elementStr.c_str());
            try
            {
                while (!elementDataStrm.fail())
                {
                    elementDataStrm >> xval >> yval >> zval;
 
                    // Need to check it here because we
                    // may not be good after the read
                    // indicating that there was nothing
                    // to read.
                    if (!elementDataStrm.fail())
                    {
                        PointGeomSharedPtr vert(
                            MemoryManager<PointGeom>::AllocateSharedPtr(
                                m_meshDimension, faceindx, xval, yval, zval));
                        curve->m_points.push_back(vert);
                    }
                }
            }
            catch (...)
            {
                NEKERROR(ErrorUtil::efatal,
                         (std::string("Unable to read curve data for FACE: ") +
                          elementStr)
                             .c_str());
            }
            m_curvedFaces[faceid] = curve;
 
            facelement = facelement->NextSiblingElement("F");
        }
    }
}

References ASSERTL0, Nektar::LibUtilities::Interpreter::DefineFunction(), Nektar::LibUtilities::Interpreter::Evaluate(), Nektar::StdRegions::find(), Nektar::LibUtilities::kPointsTypeStr, NEKERROR, and Nektar::LibUtilities::SIZE_PointsType.

v_ReadEdges()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadEdges ( )

protectedvirtual

Look for elements in ELEMENT block.

All elements are of the form: "<E ID="#"> ... </E>", with ? being the element type. Read the ID field first.

Since all edge data is one big text block, we need to accumulate all TINYXML_TEXT data and then parse it. This approach effectively skips all comments or other node types since we only care about the edge list. We cannot handle missing edge numbers as we could with missing element numbers due to the text block format.

Now parse out the edges, three fields at a time.

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 1039 of file MeshGraphXml.cpp.

{
    CurveMap::iterator it;
 
    /// Look for elements in ELEMENT block.
    TiXmlElement *field = m_xmlGeom->FirstChildElement("EDGE");
 
    ASSERTL0(field, "Unable to find EDGE tag in file.");
 
    /// All elements are of the form: "<E ID="#"> ... </E>", with
    /// ? being the element type.
    /// Read the ID field first.
    TiXmlElement *edge = field->FirstChildElement("E");
 
    /// Since all edge data is one big text block, we need to accumulate
    /// all TINYXML_TEXT data and then parse it.  This approach effectively
    /// skips
    /// all comments or other node types since we only care about the
    /// edge list.  We cannot handle missing edge numbers as we could
    /// with missing element numbers due to the text block format.
    std::string edgeStr;
    int indx;
 
    while (edge)
    {
        int err = edge->QueryIntAttribute("ID", &indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
 
        TiXmlNode *child = edge->FirstChild();
        edgeStr.clear();
        if (child->Type() == TiXmlNode::TINYXML_TEXT)
        {
            edgeStr += child->ToText()->ValueStr();
        }
 
        /// Now parse out the edges, three fields at a time.
        int vertex1, vertex2;
        std::istringstream edgeDataStrm(edgeStr.c_str());
 
        try
        {
            while (!edgeDataStrm.fail())
            {
                edgeDataStrm >> vertex1 >> vertex2;
 
                // Must check after the read because we
                // may be at the end and not know it.  If
                // we are at the end we will add a
                // duplicate of the last entry if we don't
                // check here.
                if (!edgeDataStrm.fail())
                {
                    PointGeomSharedPtr vertices[2] = {GetVertex(vertex1),
                                                      GetVertex(vertex2)};
                    SegGeomSharedPtr edge;
                    it = m_curvedEdges.find(indx);
 
                    if (it == m_curvedEdges.end())
                    {
                        edge = MemoryManager<SegGeom>::AllocateSharedPtr(
                            indx, m_spaceDimension, vertices);
                    }
                    else
                    {
                        edge = MemoryManager<SegGeom>::AllocateSharedPtr(
                            indx, m_spaceDimension, vertices, it->second);
                    }
 
                    m_segGeoms[indx] = edge;
                }
            }
        }
        catch (...)
        {
            NEKERROR(
                ErrorUtil::efatal,
                (std::string("Unable to read edge data: ") + edgeStr).c_str());
        }
 
        edge = edge->NextSiblingElement("E");
    }
}

References ASSERTL0, and NEKERROR.

v_ReadElements1D()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadElements1D ( )

protectedvirtual

Look for elements in ELEMENT block.

All elements are of the form: "<S ID = n> ... </S>", with ? being the element type.

Parse out the element components corresponding to type of element. Read two vertex numbers

Keep looking for additional segments

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 1283 of file MeshGraphXml.cpp.

{
    TiXmlElement *field = NULL;
 
    /// Look for elements in ELEMENT block.
    field = m_xmlGeom->FirstChildElement("ELEMENT");
 
    ASSERTL0(field, "Unable to find ELEMENT tag in file.");
 
    /// All elements are of the form: "<S ID = n> ... </S>", with
    /// ? being the element type.
 
    TiXmlElement *segment = field->FirstChildElement("S");
    CurveMap::iterator it;
 
    while (segment)
    {
        int indx;
        int err = segment->QueryIntAttribute("ID", &indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
 
        TiXmlNode *elementChild = segment->FirstChild();
        while (elementChild && elementChild->Type() != TiXmlNode::TINYXML_TEXT)
        {
            elementChild = elementChild->NextSibling();
        }
 
        ASSERTL0(elementChild, "Unable to read element description body.");
        std::string elementStr = elementChild->ToText()->ValueStr();
 
        /// Parse out the element components corresponding to type of
        /// element.
        /// Read two vertex numbers
        int vertex1, vertex2;
        std::istringstream elementDataStrm(elementStr.c_str());
 
        try
        {
            elementDataStrm >> vertex1;
            elementDataStrm >> vertex2;
 
            ASSERTL0(!elementDataStrm.fail(),
                     (std::string("Unable to read element data for SEGMENT: ") +
                      elementStr)
                         .c_str());
 
            PointGeomSharedPtr vertices[2] = {GetVertex(vertex1),
                                              GetVertex(vertex2)};
            SegGeomSharedPtr seg;
            it = m_curvedEdges.find(indx);
 
            if (it == m_curvedEdges.end())
            {
                seg = MemoryManager<SegGeom>::AllocateSharedPtr(
                    indx, m_spaceDimension, vertices);
                seg->SetGlobalID(indx); // Set global mesh id
            }
            else
            {
                seg = MemoryManager<SegGeom>::AllocateSharedPtr(
                    indx, m_spaceDimension, vertices, it->second);
                seg->SetGlobalID(indx); // Set global mesh id
            }
            seg->SetGlobalID(indx);
            m_segGeoms[indx] = seg;
        }
        catch (...)
        {
            NEKERROR(ErrorUtil::efatal,
                     (std::string("Unable to read element data for segment: ") +
                      elementStr)
                         .c_str());
        }
        /// Keep looking for additional segments
        segment = segment->NextSiblingElement("S");
    }
}

References ASSERTL0, and NEKERROR.

v_ReadElements2D()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadElements2D ( )

protectedvirtual

Look for elements in ELEMENT block.

All elements are of the form: "<? ID="#"> ... </?>", with ? being the element type.

Read id attribute.

Read text element description.

Parse out the element components corresponding to type of element.

Create a TriGeom to hold the new definition.

Create a QuadGeom to hold the new definition.

Keep looking

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 1361 of file MeshGraphXml.cpp.

{
    /// Look for elements in ELEMENT block.
    TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
 
    ASSERTL0(field, "Unable to find ELEMENT tag in file.");
 
    // Set up curve map for curved elements on an embedded manifold.
    CurveMap::iterator it;
 
    /// All elements are of the form: "<? ID="#"> ... </?>", with
    /// ? being the element type.
 
    TiXmlElement *element = field->FirstChildElement();
 
    while (element)
    {
        std::string elementType(element->ValueStr());
 
        ASSERTL0(
            elementType == "Q" || elementType == "T",
            (std::string("Unknown 2D element type: ") + elementType).c_str());
 
        /// Read id attribute.
        int indx;
        int err = element->QueryIntAttribute("ID", &indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
 
        it = m_curvedFaces.find(indx);
 
        /// Read text element description.
        TiXmlNode *elementChild = element->FirstChild();
        std::string elementStr;
        while (elementChild)
        {
            if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
            {
                elementStr += elementChild->ToText()->ValueStr();
            }
            elementChild = elementChild->NextSibling();
        }
 
        ASSERTL0(!elementStr.empty(),
                 "Unable to read element description body.");
 
        /// Parse out the element components corresponding to type of
        /// element.
        if (elementType == "T")
        {
            // Read three edge numbers
            int edge1, edge2, edge3;
            std::istringstream elementDataStrm(elementStr.c_str());
 
            try
            {
                elementDataStrm >> edge1;
                elementDataStrm >> edge2;
                elementDataStrm >> edge3;
 
                ASSERTL0(
                    !elementDataStrm.fail(),
                    (std::string("Unable to read element data for TRIANGLE: ") +
                     elementStr)
                        .c_str());
 
                /// Create a TriGeom to hold the new definition.
                SegGeomSharedPtr edges[TriGeom::kNedges] = {
                    GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3)};
 
                TriGeomSharedPtr trigeom;
                if (it == m_curvedFaces.end())
                {
                    trigeom =
                        MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges);
                }
                else
                {
                    trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(
                        indx, edges, it->second);
                }
                trigeom->SetGlobalID(indx);
 
                m_triGeoms[indx] = trigeom;
            }
            catch (...)
            {
                NEKERROR(
                    ErrorUtil::efatal,
                    (std::string("Unable to read element data for TRIANGLE: ") +
                     elementStr)
                        .c_str());
            }
        }
        else if (elementType == "Q")
        {
            // Read four edge numbers
            int edge1, edge2, edge3, edge4;
            std::istringstream elementDataStrm(elementStr.c_str());
 
            try
            {
                elementDataStrm >> edge1;
                elementDataStrm >> edge2;
                elementDataStrm >> edge3;
                elementDataStrm >> edge4;
 
                ASSERTL0(
                    !elementDataStrm.fail(),
                    (std::string("Unable to read element data for QUAD: ") +
                     elementStr)
                        .c_str());
 
                /// Create a QuadGeom to hold the new definition.
                SegGeomSharedPtr edges[QuadGeom::kNedges] = {
                    GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3),
                    GetSegGeom(edge4)};
 
                QuadGeomSharedPtr quadgeom;
                if (it == m_curvedFaces.end())
                {
                    quadgeom =
                        MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges);
                }
                else
                {
                    quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(
                        indx, edges, it->second);
                }
                quadgeom->SetGlobalID(indx);
 
                m_quadGeoms[indx] = quadgeom;
            }
            catch (...)
            {
                NEKERROR(
                    ErrorUtil::efatal,
                    (std::string("Unable to read element data for QUAD: ") +
                     elementStr)
                        .c_str());
            }
        }
        /// Keep looking
        element = element->NextSiblingElement();
    }
}

References ASSERTL0, and NEKERROR.

v_ReadElements3D()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadElements3D ( )

protectedvirtual

Look for elements in ELEMENT block.

All elements are of the form: "<? ID="#"> ... </?>", with ? being the element type.

Read id attribute.

Read text element description.

Parse out the element components corresponding to type of element.

Create arrays for the tri and quad faces.

Fill the arrays and make sure there aren't too many faces.

Make sure all of the face indicies could be read, and that there weren't too few.

Create arrays for the tri and quad faces.

Fill the arrays and make sure there aren't too many faces.

Make sure all of the face indicies could be read, and that there weren't too few.

Create arrays for the tri and quad faces.

Fill the arrays and make sure there aren't too many faces.

Make sure all of the face indicies could be read, and that there weren't too few.

Create arrays for the tri and quad faces.

Fill the arrays and make sure there aren't too many faces.

Make sure all of the face indicies could be read, and that there weren't too few.

Keep looking

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 1507 of file MeshGraphXml.cpp.

{
    /// Look for elements in ELEMENT block.
    TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
 
    ASSERTL0(field, "Unable to find ELEMENT tag in file.");
 
    /// All elements are of the form: "<? ID="#"> ... </?>", with
    /// ? being the element type.
 
    TiXmlElement *element = field->FirstChildElement();
 
    while (element)
    {
        std::string elementType(element->ValueStr());
 
        // A - tet, P - pyramid, R - prism, H - hex
        ASSERTL0(
            elementType == "A" || elementType == "P" || elementType == "R" ||
                elementType == "H",
            (std::string("Unknown 3D element type: ") + elementType).c_str());
 
        /// Read id attribute.
        int indx;
        int err = element->QueryIntAttribute("ID", &indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
 
        /// Read text element description.
        TiXmlNode *elementChild = element->FirstChild();
        std::string elementStr;
        while (elementChild)
        {
            if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
            {
                elementStr += elementChild->ToText()->ValueStr();
            }
            elementChild = elementChild->NextSibling();
        }
 
        ASSERTL0(!elementStr.empty(),
                 "Unable to read element description body.");
 
        std::istringstream elementDataStrm(elementStr.c_str());
 
        /// Parse out the element components corresponding to type of
        /// element.
 
        // Tetrahedral
        if (elementType == "A")
        {
            try
            {
                /// Create arrays for the tri and quad faces.
                const int kNfaces  = TetGeom::kNfaces;
                const int kNtfaces = TetGeom::kNtfaces;
                const int kNqfaces = TetGeom::kNqfaces;
                TriGeomSharedPtr tfaces[kNtfaces];
                int Ntfaces = 0;
                int Nqfaces = 0;
 
                /// Fill the arrays and make sure there aren't too many
                /// faces.
                std::stringstream errorstring;
                errorstring << "Element " << indx << " must have " << kNtfaces
                            << " triangle face(s), and " << kNqfaces
                            << " quadrilateral face(s).";
                for (int i = 0; i < kNfaces; i++)
                {
                    int faceID;
                    elementDataStrm >> faceID;
                    Geometry2DSharedPtr face = GetGeometry2D(faceID);
                    if (face == Geometry2DSharedPtr() ||
                        (face->GetShapeType() != LibUtilities::eTriangle &&
                         face->GetShapeType() != LibUtilities::eQuadrilateral))
                    {
                        std::stringstream errorstring;
                        errorstring << "Element " << indx
                                    << " has invalid face: " << faceID;
                        ASSERTL0(false, errorstring.str().c_str());
                    }
                    else if (face->GetShapeType() == LibUtilities::eTriangle)
                    {
                        ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
                        tfaces[Ntfaces++] = static_pointer_cast<TriGeom>(face);
                    }
                    else if (face->GetShapeType() ==
                             LibUtilities::eQuadrilateral)
                    {
                        ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
                    }
                }
 
                /// Make sure all of the face indicies could be read, and
                /// that there weren't too few.
                ASSERTL0(!elementDataStrm.fail(),
                         (std::string(
                              "Unable to read element data for TETRAHEDRON: ") +
                          elementStr)
                             .c_str());
                ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
                ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
 
                TetGeomSharedPtr tetgeom(
                    MemoryManager<TetGeom>::AllocateSharedPtr(indx, tfaces));
 
                m_tetGeoms[indx] = tetgeom;
                PopulateFaceToElMap(tetgeom, kNfaces);
            }
            catch (...)
            {
                NEKERROR(ErrorUtil::efatal,
                         (std::string(
                              "Unable to read element data for TETRAHEDRON: ") +
                          elementStr)
                             .c_str());
            }
        }
        // Pyramid
        else if (elementType == "P")
        {
            try
            {
                /// Create arrays for the tri and quad faces.
                const int kNfaces  = PyrGeom::kNfaces;
                const int kNtfaces = PyrGeom::kNtfaces;
                const int kNqfaces = PyrGeom::kNqfaces;
                Geometry2DSharedPtr faces[kNfaces];
                int Nfaces  = 0;
                int Ntfaces = 0;
                int Nqfaces = 0;
 
                /// Fill the arrays and make sure there aren't too many
                /// faces.
                std::stringstream errorstring;
                errorstring << "Element " << indx << " must have " << kNtfaces
                            << " triangle face(s), and " << kNqfaces
                            << " quadrilateral face(s).";
                for (int i = 0; i < kNfaces; i++)
                {
                    int faceID;
                    elementDataStrm >> faceID;
                    Geometry2DSharedPtr face = GetGeometry2D(faceID);
                    if (face == Geometry2DSharedPtr() ||
                        (face->GetShapeType() != LibUtilities::eTriangle &&
                         face->GetShapeType() != LibUtilities::eQuadrilateral))
                    {
                        std::stringstream errorstring;
                        errorstring << "Element " << indx
                                    << " has invalid face: " << faceID;
                        ASSERTL0(false, errorstring.str().c_str());
                    }
                    else if (face->GetShapeType() == LibUtilities::eTriangle)
                    {
                        ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
                        faces[Nfaces++] = static_pointer_cast<TriGeom>(face);
                        Ntfaces++;
                    }
                    else if (face->GetShapeType() ==
                             LibUtilities::eQuadrilateral)
                    {
                        ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
                        faces[Nfaces++] = static_pointer_cast<QuadGeom>(face);
                        Nqfaces++;
                    }
                }
 
                /// Make sure all of the face indicies could be read, and
                /// that there weren't too few.
                ASSERTL0(
                    !elementDataStrm.fail(),
                    (std::string("Unable to read element data for PYRAMID: ") +
                     elementStr)
                        .c_str());
                ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
                ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
 
                PyrGeomSharedPtr pyrgeom(
                    MemoryManager<PyrGeom>::AllocateSharedPtr(indx, faces));
 
                m_pyrGeoms[indx] = pyrgeom;
                PopulateFaceToElMap(pyrgeom, kNfaces);
            }
            catch (...)
            {
                NEKERROR(
                    ErrorUtil::efatal,
                    (std::string("Unable to read element data for PYRAMID: ") +
                     elementStr)
                        .c_str());
            }
        }
        // Prism
        else if (elementType == "R")
        {
            try
            {
                /// Create arrays for the tri and quad faces.
                const int kNfaces  = PrismGeom::kNfaces;
                const int kNtfaces = PrismGeom::kNtfaces;
                const int kNqfaces = PrismGeom::kNqfaces;
                Geometry2DSharedPtr faces[kNfaces];
                int Ntfaces = 0;
                int Nqfaces = 0;
                int Nfaces  = 0;
 
                /// Fill the arrays and make sure there aren't too many
                /// faces.
                std::stringstream errorstring;
                errorstring << "Element " << indx << " must have " << kNtfaces
                            << " triangle face(s), and " << kNqfaces
                            << " quadrilateral face(s).";
 
                for (int i = 0; i < kNfaces; i++)
                {
                    int faceID;
                    elementDataStrm >> faceID;
                    Geometry2DSharedPtr face = GetGeometry2D(faceID);
                    if (face == Geometry2DSharedPtr() ||
                        (face->GetShapeType() != LibUtilities::eTriangle &&
                         face->GetShapeType() != LibUtilities::eQuadrilateral))
                    {
                        std::stringstream errorstring;
                        errorstring << "Element " << indx
                                    << " has invalid face: " << faceID;
                        ASSERTL0(false, errorstring.str().c_str());
                    }
                    else if (face->GetShapeType() == LibUtilities::eTriangle)
                    {
                        ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
                        faces[Nfaces++] =
                            std::static_pointer_cast<TriGeom>(face);
                        Ntfaces++;
                    }
                    else if (face->GetShapeType() ==
                             LibUtilities::eQuadrilateral)
                    {
                        ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
                        faces[Nfaces++] =
                            std::static_pointer_cast<QuadGeom>(face);
                        Nqfaces++;
                    }
                }
 
                /// Make sure all of the face indicies could be read, and
                /// that there weren't too few.
                ASSERTL0(
                    !elementDataStrm.fail(),
                    (std::string("Unable to read element data for PRISM: ") +
                     elementStr)
                        .c_str());
                ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
                ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
 
                PrismGeomSharedPtr prismgeom(
                    MemoryManager<PrismGeom>::AllocateSharedPtr(indx, faces));
 
                m_prismGeoms[indx] = prismgeom;
                PopulateFaceToElMap(prismgeom, kNfaces);
            }
            catch (...)
            {
                NEKERROR(
                    ErrorUtil::efatal,
                    (std::string("Unable to read element data for PRISM: ") +
                     elementStr)
                        .c_str());
            }
        }
        // Hexahedral
        else if (elementType == "H")
        {
            try
            {
                /// Create arrays for the tri and quad faces.
                const int kNfaces  = HexGeom::kNfaces;
                const int kNtfaces = HexGeom::kNtfaces;
                const int kNqfaces = HexGeom::kNqfaces;
                // TriGeomSharedPtr tfaces[kNtfaces];
                QuadGeomSharedPtr qfaces[kNqfaces];
                int Ntfaces = 0;
                int Nqfaces = 0;
 
                /// Fill the arrays and make sure there aren't too many
                /// faces.
                std::stringstream errorstring;
                errorstring << "Element " << indx << " must have " << kNtfaces
                            << " triangle face(s), and " << kNqfaces
                            << " quadrilateral face(s).";
                for (int i = 0; i < kNfaces; i++)
                {
                    int faceID;
                    elementDataStrm >> faceID;
                    Geometry2DSharedPtr face = GetGeometry2D(faceID);
                    if (face == Geometry2DSharedPtr() ||
                        (face->GetShapeType() != LibUtilities::eTriangle &&
                         face->GetShapeType() != LibUtilities::eQuadrilateral))
                    {
                        std::stringstream errorstring;
                        errorstring << "Element " << indx
                                    << " has invalid face: " << faceID;
                        ASSERTL0(false, errorstring.str().c_str());
                    }
                    else if (face->GetShapeType() == LibUtilities::eTriangle)
                    {
                        ASSERTL0(Ntfaces < kNtfaces, errorstring.str().c_str());
                        // tfaces[Ntfaces++] =
                        // boost::static_pointer_cast<TriGeom>(face);
                    }
                    else if (face->GetShapeType() ==
                             LibUtilities::eQuadrilateral)
                    {
                        ASSERTL0(Nqfaces < kNqfaces, errorstring.str().c_str());
                        qfaces[Nqfaces++] =
                            std::static_pointer_cast<QuadGeom>(face);
                    }
                }
 
                /// Make sure all of the face indicies could be read, and
                /// that there weren't too few.
                ASSERTL0(!elementDataStrm.fail(),
                         (std::string(
                              "Unable to read element data for HEXAHEDRAL: ") +
                          elementStr)
                             .c_str());
                ASSERTL0(Ntfaces == kNtfaces, errorstring.str().c_str());
                ASSERTL0(Nqfaces == kNqfaces, errorstring.str().c_str());
 
                HexGeomSharedPtr hexgeom(
                    MemoryManager<HexGeom>::AllocateSharedPtr(indx, qfaces));
 
                m_hexGeoms[indx] = hexgeom;
                PopulateFaceToElMap(hexgeom, kNfaces);
            }
            catch (...)
            {
                NEKERROR(ErrorUtil::efatal,
                         (std::string(
                              "Unable to read element data for HEXAHEDRAL: ") +
                          elementStr)
                             .c_str());
            }
        }
        /// Keep looking
        element = element->NextSiblingElement();
    }
}

References ASSERTL0, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, and NEKERROR.

v_ReadFaces()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadFaces ( )

protectedvirtual

Look for elements in FACE block.

All faces are of the form: "<? ID="#"> ... </?>", with ? being an element type (either Q or T). They might be in compressed format and so then need upacking.

Read id attribute.

See if this face has curves.

Read text element description.

Parse out the element components corresponding to type of element.

Create a TriGeom to hold the new definition.

Create a QuadGeom to hold the new definition.

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 1122 of file MeshGraphXml.cpp.

{
    /// Look for elements in FACE block.
    TiXmlElement *field = m_xmlGeom->FirstChildElement("FACE");
 
    ASSERTL0(field, "Unable to find FACE tag in file.");
 
    /// All faces are of the form: "<? ID="#"> ... </?>", with
    /// ? being an element type (either Q or T).
    /// They might be in compressed format and so then need upacking.
 
    TiXmlElement *element = field->FirstChildElement();
    CurveMap::iterator it;
 
    while (element)
    {
        std::string elementType(element->ValueStr());
 
        ASSERTL0(elementType == "Q" || elementType == "T",
                 (std::string("Unknown 3D face type: ") + elementType).c_str());
 
        /// Read id attribute.
        int indx;
        int err = element->QueryIntAttribute("ID", &indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read face attribute ID.");
 
        /// See if this face has curves.
        it = m_curvedFaces.find(indx);
 
        /// Read text element description.
        TiXmlNode *elementChild = element->FirstChild();
        std::string elementStr;
        while (elementChild)
        {
            if (elementChild->Type() == TiXmlNode::TINYXML_TEXT)
            {
                elementStr += elementChild->ToText()->ValueStr();
            }
            elementChild = elementChild->NextSibling();
        }
 
        ASSERTL0(!elementStr.empty(), "Unable to read face description body.");
 
        /// Parse out the element components corresponding to type of
        /// element.
        if (elementType == "T")
        {
            // Read three edge numbers
            int edge1, edge2, edge3;
            std::istringstream elementDataStrm(elementStr.c_str());
 
            try
            {
                elementDataStrm >> edge1;
                elementDataStrm >> edge2;
                elementDataStrm >> edge3;
 
                ASSERTL0(
                    !elementDataStrm.fail(),
                    (std::string("Unable to read face data for TRIANGLE: ") +
                     elementStr)
                        .c_str());
 
                /// Create a TriGeom to hold the new definition.
                SegGeomSharedPtr edges[TriGeom::kNedges] = {
                    GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3)};
 
                TriGeomSharedPtr trigeom;
 
                if (it == m_curvedFaces.end())
                {
                    trigeom =
                        MemoryManager<TriGeom>::AllocateSharedPtr(indx, edges);
                }
                else
                {
                    trigeom = MemoryManager<TriGeom>::AllocateSharedPtr(
                        indx, edges, it->second);
                }
 
                trigeom->SetGlobalID(indx);
 
                m_triGeoms[indx] = trigeom;
            }
            catch (...)
            {
                NEKERROR(
                    ErrorUtil::efatal,
                    (std::string("Unable to read face data for TRIANGLE: ") +
                     elementStr)
                        .c_str());
            }
        }
        else if (elementType == "Q")
        {
            // Read four edge numbers
            int edge1, edge2, edge3, edge4;
            std::istringstream elementDataStrm(elementStr.c_str());
 
            try
            {
                elementDataStrm >> edge1;
                elementDataStrm >> edge2;
                elementDataStrm >> edge3;
                elementDataStrm >> edge4;
 
                ASSERTL0(!elementDataStrm.fail(),
                         (std::string("Unable to read face data for QUAD: ") +
                          elementStr)
                             .c_str());
 
                /// Create a QuadGeom to hold the new definition.
                SegGeomSharedPtr edges[QuadGeom::kNedges] = {
                    GetSegGeom(edge1), GetSegGeom(edge2), GetSegGeom(edge3),
                    GetSegGeom(edge4)};
 
                QuadGeomSharedPtr quadgeom;
 
                if (it == m_curvedFaces.end())
                {
                    quadgeom =
                        MemoryManager<QuadGeom>::AllocateSharedPtr(indx, edges);
                }
                else
                {
                    quadgeom = MemoryManager<QuadGeom>::AllocateSharedPtr(
                        indx, edges, it->second);
                }
                quadgeom->SetGlobalID(indx);
 
                m_quadGeoms[indx] = quadgeom;
            }
            catch (...)
            {
                NEKERROR(ErrorUtil::efatal,
                         (std::string("Unable to read face data for QUAD: ") +
                          elementStr)
                             .c_str());
            }
        }
        // Keep looking
        element = element->NextSiblingElement();
    }
}

References ASSERTL0, and NEKERROR.

v_ReadGeometry()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadGeometry ( LibUtilities::DomainRangeShPtr  rng, bool  fillGraph  )

overrideprotectedvirtual

Error value returned by TinyXML.

Implements Nektar::SpatialDomains::MeshGraph.

Definition at line 376 of file MeshGraphXml.cpp.

{
    // Reset member variables.
    m_vertSet.clear();
    m_curvedEdges.clear();
    m_curvedFaces.clear();
    m_segGeoms.clear();
    m_triGeoms.clear();
    m_quadGeoms.clear();
    m_tetGeoms.clear();
    m_pyrGeoms.clear();
    m_prismGeoms.clear();
    m_hexGeoms.clear();
    m_meshComposites.clear();
    m_compositesLabels.clear();
    m_domain.clear();
    m_expansionMapShPtrMap.clear();
    m_geomInfo.clear();
    m_faceToElMap.clear();
 
    m_domainRange = rng;
    m_xmlGeom     = m_session->GetElement("NEKTAR/GEOMETRY");
 
    int err; /// Error value returned by TinyXML.
 
    TiXmlAttribute *attr = m_xmlGeom->FirstAttribute();
 
    // Initialize the mesh and space dimensions to 3 dimensions.
    // We want to do this each time we read a file, so it should
    // be done here and not just during class initialization.
    m_meshPartitioned = false;
    m_meshDimension   = 3;
    m_spaceDimension  = 3;
 
    while (attr)
    {
        std::string attrName(attr->Name());
        if (attrName == "DIM")
        {
            err = attr->QueryIntValue(&m_meshDimension);
            ASSERTL0(err == TIXML_SUCCESS, "Unable to read mesh dimension.");
        }
        else if (attrName == "SPACE")
        {
            err = attr->QueryIntValue(&m_spaceDimension);
            ASSERTL0(err == TIXML_SUCCESS, "Unable to read space dimension.");
        }
        else if (attrName == "PARTITION")
        {
            err = attr->QueryIntValue(&m_partition);
            ASSERTL0(err == TIXML_SUCCESS, "Unable to read partition.");
            m_meshPartitioned = true;
        }
        else
        {
            std::string errstr("Unknown attribute: ");
            errstr += attrName;
            ASSERTL0(false, errstr.c_str());
        }
 
        // Get the next attribute.
        attr = attr->Next();
    }
 
    ASSERTL0(m_meshDimension <= m_spaceDimension,
             "Mesh dimension greater than space dimension");
 
    v_ReadVertices();
    v_ReadCurves();
    if (m_meshDimension >= 2)
    {
        v_ReadEdges();
        if (m_meshDimension == 3)
        {
            v_ReadFaces();
        }
    }
    ReadElements();
    ReadComposites();
    ReadDomain();
 
    if (fillGraph)
    {
        MeshGraph::FillGraph();
    }
}

References ASSERTL0.

v_ReadVertices()

void Nektar::SpatialDomains::MeshGraphXml::v_ReadVertices ( )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 464 of file MeshGraphXml.cpp.

{
    // Now read the vertices
    TiXmlElement *element = m_xmlGeom->FirstChildElement("VERTEX");
    ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
 
    NekDouble xscale, yscale, zscale;
 
    // check to see if any scaling parameters are in
    // attributes and determine these values
    LibUtilities::Interpreter expEvaluator;
    const char *xscal = element->Attribute("XSCALE");
    if (!xscal)
    {
        xscale = 1.0;
    }
    else
    {
        std::string xscalstr = xscal;
        int expr_id          = expEvaluator.DefineFunction("", xscalstr);
        xscale               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *yscal = element->Attribute("YSCALE");
    if (!yscal)
    {
        yscale = 1.0;
    }
    else
    {
        std::string yscalstr = yscal;
        int expr_id          = expEvaluator.DefineFunction("", yscalstr);
        yscale               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *zscal = element->Attribute("ZSCALE");
    if (!zscal)
    {
        zscale = 1.0;
    }
    else
    {
        std::string zscalstr = zscal;
        int expr_id          = expEvaluator.DefineFunction("", zscalstr);
        zscale               = expEvaluator.Evaluate(expr_id);
    }
 
    NekDouble xmove, ymove, zmove;
 
    // check to see if any moving parameters are in
    // attributes and determine these values
 
    const char *xmov = element->Attribute("XMOVE");
    if (!xmov)
    {
        xmove = 0.0;
    }
    else
    {
        std::string xmovstr = xmov;
        int expr_id         = expEvaluator.DefineFunction("", xmovstr);
        xmove               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *ymov = element->Attribute("YMOVE");
    if (!ymov)
    {
        ymove = 0.0;
    }
    else
    {
        std::string ymovstr = ymov;
        int expr_id         = expEvaluator.DefineFunction("", ymovstr);
        ymove               = expEvaluator.Evaluate(expr_id);
    }
 
    const char *zmov = element->Attribute("ZMOVE");
    if (!zmov)
    {
        zmove = 0.0;
    }
    else
    {
        std::string zmovstr = zmov;
        int expr_id         = expEvaluator.DefineFunction("", zmovstr);
        zmove               = expEvaluator.Evaluate(expr_id);
    }
 
    TiXmlElement *vertex = element->FirstChildElement("V");
 
    int indx;
    int nextVertexNumber = -1;
 
    while (vertex)
    {
        nextVertexNumber++;
 
        TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
        std::string attrName(vertexAttr->Name());
 
        ASSERTL0(attrName == "ID",
                 (std::string("Unknown attribute name: ") + attrName).c_str());
 
        int err = vertexAttr->QueryIntValue(&indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
 
        // Now read body of vertex
        std::string vertexBodyStr;
 
        TiXmlNode *vertexBody = vertex->FirstChild();
 
        while (vertexBody)
        {
            // Accumulate all non-comment body data.
            if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
            {
                vertexBodyStr += vertexBody->ToText()->Value();
                vertexBodyStr += " ";
            }
 
            vertexBody = vertexBody->NextSibling();
        }
 
        ASSERTL0(!vertexBodyStr.empty(),
                 "Vertex definitions must contain vertex data.");
 
        // Get vertex data from the data string.
        NekDouble xval, yval, zval;
        std::istringstream vertexDataStrm(vertexBodyStr.c_str());
 
        try
        {
            while (!vertexDataStrm.fail())
            {
                vertexDataStrm >> xval >> yval >> zval;
 
                xval = xval * xscale + xmove;
                yval = yval * yscale + ymove;
                zval = zval * zscale + zmove;
 
                // Need to check it here because we may not be
                // good after the read indicating that there
                // was nothing to read.
                if (!vertexDataStrm.fail())
                {
                    PointGeomSharedPtr vert(
                        MemoryManager<PointGeom>::AllocateSharedPtr(
                            m_spaceDimension, indx, xval, yval, zval));
                    m_vertSet[indx] = vert;
                }
            }
        }
        catch (...)
        {
            ASSERTL0(false, "Unable to read VERTEX data.");
        }
 
        vertex = vertex->NextSiblingElement("V");
    }
}

References ASSERTL0, Nektar::LibUtilities::Interpreter::DefineFunction(), and Nektar::LibUtilities::Interpreter::Evaluate().

v_WriteCurves()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteCurves ( TiXmlElement *  geomTag, CurveMap &  edges, CurveMap &  faces  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2603 of file MeshGraphXml.cpp.

{
    TiXmlElement *curveTag = new TiXmlElement("CURVED");
    CurveMap::iterator curveIt;
    int curveId = 0;
 
    for (curveIt = edges.begin(); curveIt != edges.end(); ++curveIt)
    {
        CurveSharedPtr curve = curveIt->second;
        TiXmlElement *c      = new TiXmlElement("E");
        stringstream s;
        s.precision(8);
 
        for (int j = 0; j < curve->m_points.size(); ++j)
        {
            SpatialDomains::PointGeomSharedPtr p = curve->m_points[j];
            s << scientific << (*p)(0) << " " << (*p)(1) << " " << (*p)(2)
              << "   ";
        }
 
        c->SetAttribute("ID", curveId++);
        c->SetAttribute("EDGEID", curve->m_curveID);
        c->SetAttribute("NUMPOINTS", curve->m_points.size());
        c->SetAttribute("TYPE", LibUtilities::kPointsTypeStr[curve->m_ptype]);
        c->LinkEndChild(new TiXmlText(s.str()));
        curveTag->LinkEndChild(c);
    }
 
    for (curveIt = faces.begin(); curveIt != faces.end(); ++curveIt)
    {
        CurveSharedPtr curve = curveIt->second;
        TiXmlElement *c      = new TiXmlElement("F");
        stringstream s;
        s.precision(8);
 
        for (int j = 0; j < curve->m_points.size(); ++j)
        {
            SpatialDomains::PointGeomSharedPtr p = curve->m_points[j];
            s << scientific << (*p)(0) << " " << (*p)(1) << " " << (*p)(2)
              << "   ";
        }
 
        c->SetAttribute("ID", curveId++);
        c->SetAttribute("FACEID", curve->m_curveID);
        c->SetAttribute("NUMPOINTS", curve->m_points.size());
        c->SetAttribute("TYPE", LibUtilities::kPointsTypeStr[curve->m_ptype]);
        c->LinkEndChild(new TiXmlText(s.str()));
        curveTag->LinkEndChild(c);
    }
 
    geomTag->LinkEndChild(curveTag);
}

References Nektar::LibUtilities::kPointsTypeStr, and CellMLToNektar.cellml_metadata::p.

v_WriteEdges()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteEdges ( TiXmlElement *  geomTag, SegGeomMap &  edges  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2480 of file MeshGraphXml.cpp.

{
    TiXmlElement *edgeTag =
        new TiXmlElement(m_meshDimension == 1 ? "ELEMENT" : "EDGE");
    string tag = m_meshDimension == 1 ? "S" : "E";
 
    for (auto &i : edges)
    {
        stringstream s;
        SegGeomSharedPtr seg = i.second;
        s << seg->GetVid(0) << " " << seg->GetVid(1);
        TiXmlElement *e = new TiXmlElement(tag);
        e->SetAttribute("ID", i.first);
        e->LinkEndChild(new TiXmlText(s.str()));
        edgeTag->LinkEndChild(e);
    }
 
    geomTag->LinkEndChild(edgeTag);
}

v_WriteGeometry()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteGeometry ( std::string &  outfilename, bool  defaultExp = false, const LibUtilities::FieldMetaDataMap &  metadata = LibUtilities::NullFieldMetaDataMap  )

overrideprotectedvirtual

Write out an XML file containing the GEOMETRY block representing this MeshGraph instance inside a NEKTAR tag.

Implements Nektar::SpatialDomains::MeshGraph.

Definition at line 2737 of file MeshGraphXml.cpp.

{
    // Create empty TinyXML document.
    TiXmlDocument doc;
    TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", "");
    doc.LinkEndChild(decl);
 
    TiXmlElement *root = new TiXmlElement("NEKTAR");
    doc.LinkEndChild(root);
    TiXmlElement *geomTag = new TiXmlElement("GEOMETRY");
    root->LinkEndChild(geomTag);
 
    // Add provenance information using FieldIO library.
    LibUtilities::FieldIO::AddInfoTag(LibUtilities::XmlTagWriterSharedPtr(
                                          new LibUtilities::XmlTagWriter(root)),
                                      metadata);
 
    // Update attributes with dimensions.
    geomTag->SetAttribute("DIM", m_meshDimension);
    geomTag->SetAttribute("SPACE", m_spaceDimension);
 
    if (m_session != nullptr && !m_session->GetComm()->IsSerial())
    {
        geomTag->SetAttribute("PARTITION", m_session->GetComm()->GetRank());
    }
 
    // Clear existing elements.
    geomTag->Clear();
 
    // Write out informatio
    v_WriteVertices(geomTag, m_vertSet);
    v_WriteEdges(geomTag, m_segGeoms);
    if (m_meshDimension > 1)
    {
        TiXmlElement *faceTag =
            new TiXmlElement(m_meshDimension == 2 ? "ELEMENT" : "FACE");
 
        v_WriteTris(faceTag, m_triGeoms);
        v_WriteQuads(faceTag, m_quadGeoms);
        geomTag->LinkEndChild(faceTag);
    }
    if (m_meshDimension > 2)
    {
        TiXmlElement *elmtTag = new TiXmlElement("ELEMENT");
 
        v_WriteHexs(elmtTag, m_hexGeoms);
        v_WritePyrs(elmtTag, m_pyrGeoms);
        v_WritePrisms(elmtTag, m_prismGeoms);
        v_WriteTets(elmtTag, m_tetGeoms);
 
        geomTag->LinkEndChild(elmtTag);
    }
    v_WriteCurves(geomTag, m_curvedEdges, m_curvedFaces);
    WriteComposites(geomTag, m_meshComposites, m_compositesLabels);
    WriteDomain(geomTag, m_domain);
 
    if (defaultExp)
    {
        WriteDefaultExpansion(root);
    }
 
    // Save file.
    doc.SaveFile(outfilename);
}

v_WriteHexs()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteHexs ( TiXmlElement *  elmtTag, HexGeomMap &  hexs  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2533 of file MeshGraphXml.cpp.

{
    string tag = "H";
 
    for (auto &i : hexs)
    {
        stringstream s;
        HexGeomSharedPtr hex = i.second;
        s << hex->GetFid(0) << " " << hex->GetFid(1) << " " << hex->GetFid(2)
          << " " << hex->GetFid(3) << " " << hex->GetFid(4) << " "
          << hex->GetFid(5) << " ";
        TiXmlElement *h = new TiXmlElement(tag);
        h->SetAttribute("ID", i.first);
        h->LinkEndChild(new TiXmlText(s.str()));
        elmtTag->LinkEndChild(h);
    }
}

v_WritePrisms()

void Nektar::SpatialDomains::MeshGraphXml::v_WritePrisms ( TiXmlElement *  elmtTag, PrismGeomMap &  pris  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2551 of file MeshGraphXml.cpp.

{
    string tag = "R";
 
    for (auto &i : pris)
    {
        stringstream s;
        PrismGeomSharedPtr prism = i.second;
        s << prism->GetFid(0) << " " << prism->GetFid(1) << " "
          << prism->GetFid(2) << " " << prism->GetFid(3) << " "
          << prism->GetFid(4) << " ";
        TiXmlElement *p = new TiXmlElement(tag);
        p->SetAttribute("ID", i.first);
        p->LinkEndChild(new TiXmlText(s.str()));
        elmtTag->LinkEndChild(p);
    }
}

References CellMLToNektar.cellml_metadata::p.

v_WritePyrs()

void Nektar::SpatialDomains::MeshGraphXml::v_WritePyrs ( TiXmlElement *  elmtTag, PyrGeomMap &  pyrs  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2569 of file MeshGraphXml.cpp.

{
    string tag = "P";
 
    for (auto &i : pyrs)
    {
        stringstream s;
        PyrGeomSharedPtr pyr = i.second;
        s << pyr->GetFid(0) << " " << pyr->GetFid(1) << " " << pyr->GetFid(2)
          << " " << pyr->GetFid(3) << " " << pyr->GetFid(4) << " ";
        TiXmlElement *p = new TiXmlElement(tag);
        p->SetAttribute("ID", i.first);
        p->LinkEndChild(new TiXmlText(s.str()));
        elmtTag->LinkEndChild(p);
    }
}

References CellMLToNektar.cellml_metadata::p.

v_WriteQuads()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteQuads ( TiXmlElement *  faceTag, QuadGeomMap &  quads  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2516 of file MeshGraphXml.cpp.

{
    string tag = "Q";
 
    for (auto &i : quads)
    {
        stringstream s;
        QuadGeomSharedPtr quad = i.second;
        s << quad->GetEid(0) << " " << quad->GetEid(1) << " " << quad->GetEid(2)
          << " " << quad->GetEid(3);
        TiXmlElement *q = new TiXmlElement(tag);
        q->SetAttribute("ID", i.first);
        q->LinkEndChild(new TiXmlText(s.str()));
        faceTag->LinkEndChild(q);
    }
}

v_WriteTets()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteTets ( TiXmlElement *  elmtTag, TetGeomMap &  tets  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2586 of file MeshGraphXml.cpp.

{
    string tag = "A";
 
    for (auto &i : tets)
    {
        stringstream s;
        TetGeomSharedPtr tet = i.second;
        s << tet->GetFid(0) << " " << tet->GetFid(1) << " " << tet->GetFid(2)
          << " " << tet->GetFid(3) << " ";
        TiXmlElement *t = new TiXmlElement(tag);
        t->SetAttribute("ID", i.first);
        t->LinkEndChild(new TiXmlText(s.str()));
        elmtTag->LinkEndChild(t);
    }
}

v_WriteTris()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteTris ( TiXmlElement *  faceTag, TriGeomMap &  tris  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2500 of file MeshGraphXml.cpp.

{
    string tag = "T";
 
    for (auto &i : tris)
    {
        stringstream s;
        TriGeomSharedPtr tri = i.second;
        s << tri->GetEid(0) << " " << tri->GetEid(1) << " " << tri->GetEid(2);
        TiXmlElement *t = new TiXmlElement(tag);
        t->SetAttribute("ID", i.first);
        t->LinkEndChild(new TiXmlText(s.str()));
        faceTag->LinkEndChild(t);
    }
}

v_WriteVertices()

void Nektar::SpatialDomains::MeshGraphXml::v_WriteVertices ( TiXmlElement *  geomTag, PointGeomMap &  verts  )

protectedvirtual

Reimplemented in Nektar::SpatialDomains::MeshGraphXmlCompressed.

Definition at line 2462 of file MeshGraphXml.cpp.

{
    TiXmlElement *vertTag = new TiXmlElement("VERTEX");
 
    for (auto &i : verts)
    {
        stringstream s;
        s << scientific << setprecision(8) << (*i.second)(0) << " "
          << (*i.second)(1) << " " << (*i.second)(2);
        TiXmlElement *v = new TiXmlElement("V");
        v->SetAttribute("ID", i.second->GetGlobalID());
        v->LinkEndChild(new TiXmlText(s.str()));
        vertTag->LinkEndChild(v);
    }
 
    geomTag->LinkEndChild(vertTag);
}

WriteComposites()

void Nektar::SpatialDomains::MeshGraphXml::WriteComposites ( TiXmlElement *  geomTag, CompositeMap &  comps, std::map< int, std::string > &  compLabels  )

protected

Definition at line 2657 of file MeshGraphXml.cpp.

{
    TiXmlElement *compTag = new TiXmlElement("COMPOSITE");
 
    for (auto &cIt : comps)
    {
        if (cIt.second->m_geomVec.size() == 0)
        {
            continue;
        }
 
        TiXmlElement *c = new TiXmlElement("C");
        c->SetAttribute("ID", cIt.first);
        if (!m_compositesLabels[cIt.first].empty())
        {
            c->SetAttribute("NAME", compLabels[cIt.first]);
        }
        c->LinkEndChild(new TiXmlText(GetCompositeString(cIt.second)));
        compTag->LinkEndChild(c);
    }
 
    geomTag->LinkEndChild(compTag);
}

WriteDefaultExpansion()

void Nektar::SpatialDomains::MeshGraphXml::WriteDefaultExpansion ( TiXmlElement *  root )

protected

Definition at line 2711 of file MeshGraphXml.cpp.

{
    TiXmlElement *expTag = new TiXmlElement("EXPANSIONS");
 
    for (auto it = m_meshComposites.begin(); it != m_meshComposites.end(); it++)
    {
        if (it->second->m_geomVec[0]->GetShapeDim() == m_meshDimension)
        {
            TiXmlElement *exp = new TiXmlElement("E");
            exp->SetAttribute("COMPOSITE",
                              "C[" + boost::lexical_cast<string>(it->first) +
                                  "]");
            exp->SetAttribute("NUMMODES", 4);
            exp->SetAttribute("TYPE", "MODIFIED");
            exp->SetAttribute("FIELDS", "u");
 
            expTag->LinkEndChild(exp);
        }
    }
    root->LinkEndChild(expTag);
}

WriteDomain()

void Nektar::SpatialDomains::MeshGraphXml::WriteDomain ( TiXmlElement *  geomTag, std::map< int, CompositeMap > &  domain  )

protected

Definition at line 2682 of file MeshGraphXml.cpp.

{
    TiXmlElement *domTag = new TiXmlElement("DOMAIN");
 
    vector<unsigned int> idxList;
    for (auto &domain : domainMap)
    {
        TiXmlElement *c = new TiXmlElement("D");
        idxList.clear();
        stringstream s;
        s << " "
          << "C"
          << "[";
 
        for (const auto &elem : domain.second)
        {
            idxList.push_back(elem.first);
        }
 
        s << ParseUtils::GenerateSeqString(idxList) << "] ";
        c->SetAttribute("ID", domain.first);
        c->LinkEndChild(new TiXmlText(s.str()));
        domTag->LinkEndChild(c);
    }
 
    geomTag->LinkEndChild(domTag);
}

WriteXMLGeometry()

void Nektar::SpatialDomains::MeshGraphXml::WriteXMLGeometry	(	std::string	outname,
		std::vector< std::set< unsigned int >>	elements,
		std::vector< unsigned int >	partitions
	)

Definition at line 2804 of file MeshGraphXml.cpp.

{
    // so in theory this function is used by the mesh partitioner
    // giving instructions on how to write out a paritioned mesh.
    // the theory goes that the elements stored in the meshgraph are the
    // "whole" mesh so based on the information from the elmements list
    // we can filter the mesh entities and write some individual files
    // hopefully
 
    // this is xml so we are going to write a directory with lots of
    // xml files
    string dirname = outname + "_xml";
    boost::filesystem::path pdirname(dirname);
 
    if (!boost::filesystem::is_directory(dirname))
    {
        boost::filesystem::create_directory(dirname);
    }
 
    ASSERTL0(elements.size() == partitions.size(),
             "error in partitioned information");
 
    for (int i = 0; i < partitions.size(); i++)
    {
        TiXmlDocument doc;
        TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", "");
        doc.LinkEndChild(decl);
 
        TiXmlElement *root = doc.FirstChildElement("NEKTAR");
        TiXmlElement *geomTag;
 
        // Try to find existing NEKTAR tag.
        if (!root)
        {
            root = new TiXmlElement("NEKTAR");
            doc.LinkEndChild(root);
 
            geomTag = new TiXmlElement("GEOMETRY");
            root->LinkEndChild(geomTag);
        }
        else
        {
            // Try to find existing GEOMETRY tag.
            geomTag = root->FirstChildElement("GEOMETRY");
 
            if (!geomTag)
            {
                geomTag = new TiXmlElement("GEOMETRY");
                root->LinkEndChild(geomTag);
            }
        }
 
        geomTag->SetAttribute("DIM", m_meshDimension);
        geomTag->SetAttribute("SPACE", m_spaceDimension);
        geomTag->SetAttribute("PARTITION", partitions[i]);
 
        // Add Mesh //
        // Get the elements
        HexGeomMap localHex;
        PyrGeomMap localPyr;
        PrismGeomMap localPrism;
        TetGeomMap localTet;
        TriGeomMap localTri;
        QuadGeomMap localQuad;
        SegGeomMap localEdge;
        PointGeomMap localVert;
        CurveMap localCurveEdge;
        CurveMap localCurveFace;
 
        vector<set<unsigned int>> entityIds(4);
        entityIds[m_meshDimension] = elements[i];
 
        switch (m_meshDimension)
        {
            case 3:
            {
                for (auto &j : entityIds[3])
                {
                    GeometrySharedPtr g;
                    if (m_hexGeoms.count(j))
                    {
                        g           = m_hexGeoms[j];
                        localHex[j] = m_hexGeoms[j];
                    }
                    else if (m_pyrGeoms.count(j))
                    {
                        g           = m_pyrGeoms[j];
                        localPyr[j] = m_pyrGeoms[j];
                    }
                    else if (m_prismGeoms.count(j))
                    {
                        g             = m_prismGeoms[j];
                        localPrism[j] = m_prismGeoms[j];
                    }
                    else if (m_tetGeoms.count(j))
                    {
                        g           = m_tetGeoms[j];
                        localTet[j] = m_tetGeoms[j];
                    }
                    else
                    {
                        ASSERTL0(false, "element in partition not found");
                    }
 
                    for (int k = 0; k < g->GetNumFaces(); k++)
                    {
                        entityIds[2].insert(g->GetFid(k));
                    }
                    for (int k = 0; k < g->GetNumEdges(); k++)
                    {
                        entityIds[1].insert(g->GetEid(k));
                    }
                    for (int k = 0; k < g->GetNumVerts(); k++)
                    {
                        entityIds[0].insert(g->GetVid(k));
                    }
                }
            }
            break;
            case 2:
            {
                for (auto &j : entityIds[2])
                {
                    GeometrySharedPtr g;
                    if (m_triGeoms.count(j))
                    {
                        g           = m_triGeoms[j];
                        localTri[j] = m_triGeoms[j];
                    }
                    else if (m_quadGeoms.count(j))
                    {
                        g            = m_quadGeoms[j];
                        localQuad[j] = m_quadGeoms[j];
                    }
                    else
                    {
                        ASSERTL0(false, "element in partition not found");
                    }
 
                    for (int k = 0; k < g->GetNumEdges(); k++)
                    {
                        entityIds[1].insert(g->GetEid(k));
                    }
                    for (int k = 0; k < g->GetNumVerts(); k++)
                    {
                        entityIds[0].insert(g->GetVid(k));
                    }
                }
            }
            break;
            case 1:
            {
                for (auto &j : entityIds[1])
                {
                    GeometrySharedPtr g;
                    if (m_segGeoms.count(j))
                    {
                        g            = m_segGeoms[j];
                        localEdge[j] = m_segGeoms[j];
                    }
                    else
                    {
                        ASSERTL0(false, "element in partition not found");
                    }
 
                    for (int k = 0; k < g->GetNumVerts(); k++)
                    {
                        entityIds[0].insert(g->GetVid(k));
                    }
                }
            }
        }
 
        if (m_meshDimension > 2)
        {
            for (auto &j : entityIds[2])
            {
                if (m_triGeoms.count(j))
                {
                    localTri[j] = m_triGeoms[j];
                }
                else if (m_quadGeoms.count(j))
                {
                    localQuad[j] = m_quadGeoms[j];
                }
                else
                {
                    ASSERTL0(false, "face not found");
                }
            }
        }
 
        if (m_meshDimension > 1)
        {
            for (auto &j : entityIds[1])
            {
                if (m_segGeoms.count(j))
                {
                    localEdge[j] = m_segGeoms[j];
                }
                else
                {
                    ASSERTL0(false, "edge not found");
                }
            }
        }
 
        for (auto &j : entityIds[0])
        {
            if (m_vertSet.count(j))
            {
                localVert[j] = m_vertSet[j];
            }
            else
            {
                ASSERTL0(false, "vert not found");
            }
        }
 
        v_WriteVertices(geomTag, localVert);
        v_WriteEdges(geomTag, localEdge);
        if (m_meshDimension > 1)
        {
            TiXmlElement *faceTag =
                new TiXmlElement(m_meshDimension == 2 ? "ELEMENT" : "FACE");
 
            v_WriteTris(faceTag, localTri);
            v_WriteQuads(faceTag, localQuad);
            geomTag->LinkEndChild(faceTag);
        }
        if (m_meshDimension > 2)
        {
            TiXmlElement *elmtTag = new TiXmlElement("ELEMENT");
 
            v_WriteHexs(elmtTag, localHex);
            v_WritePyrs(elmtTag, localPyr);
            v_WritePrisms(elmtTag, localPrism);
            v_WriteTets(elmtTag, localTet);
 
            geomTag->LinkEndChild(elmtTag);
        }
 
        for (auto &j : localTri)
        {
            if (m_curvedFaces.count(j.first))
            {
                localCurveFace[j.first] = m_curvedFaces[j.first];
            }
        }
        for (auto &j : localQuad)
        {
            if (m_curvedFaces.count(j.first))
            {
                localCurveFace[j.first] = m_curvedFaces[j.first];
            }
        }
        for (auto &j : localEdge)
        {
            if (m_curvedEdges.count(j.first))
            {
                localCurveEdge[j.first] = m_curvedEdges[j.first];
            }
        }
 
        v_WriteCurves(geomTag, localCurveEdge, localCurveFace);
 
        CompositeMap localComp;
        std::map<int, std::string> localCompLabels;
 
        for (auto &j : m_meshComposites)
        {
            CompositeSharedPtr comp = CompositeSharedPtr(new Composite);
            int dim                 = j.second->m_geomVec[0]->GetShapeDim();
 
            for (int k = 0; k < j.second->m_geomVec.size(); k++)
            {
                if (entityIds[dim].count(j.second->m_geomVec[k]->GetGlobalID()))
                {
                    comp->m_geomVec.push_back(j.second->m_geomVec[k]);
                }
            }
 
            if (comp->m_geomVec.size())
            {
                localComp[j.first] = comp;
                if (!m_compositesLabels[j.first].empty())
                {
                    localCompLabels[j.first] = m_compositesLabels[j.first];
                }
            }
        }
 
        WriteComposites(geomTag, localComp, localCompLabels);
 
        map<int, CompositeMap> domain;
        for (auto &j : localComp)
        {
            for (auto &dom : m_domain)
            {
                for (auto &dIt : dom.second)
                {
                    if (j.first == dIt.first)
                    {
                        domain[dom.first][j.first] = j.second;
                        break;
                    }
                }
            }
        }
 
        WriteDomain(geomTag, domain);
 
        if (m_session->DefinesElement("NEKTAR/CONDITIONS"))
        {
            std::set<int> vBndRegionIdList;
            TiXmlElement *vConditions =
                new TiXmlElement(*m_session->GetElement("Nektar/Conditions"));
            TiXmlElement *vBndRegions =
                vConditions->FirstChildElement("BOUNDARYREGIONS");
            TiXmlElement *vBndConditions =
                vConditions->FirstChildElement("BOUNDARYCONDITIONS");
            TiXmlElement *vItem;
 
            if (vBndRegions)
            {
                TiXmlElement *vNewBndRegions =
                    new TiXmlElement("BOUNDARYREGIONS");
                vItem = vBndRegions->FirstChildElement();
                while (vItem)
                {
                    std::string vSeqStr =
                        vItem->FirstChild()->ToText()->Value();
                    std::string::size_type indxBeg =
                        vSeqStr.find_first_of('[') + 1;
                    std::string::size_type indxEnd =
                        vSeqStr.find_last_of(']') - 1;
                    vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
                    std::vector<unsigned int> vSeq;
                    ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
 
                    vector<unsigned int> idxList;
 
                    for (unsigned int i = 0; i < vSeq.size(); ++i)
                    {
                        if (localComp.find(vSeq[i]) != localComp.end())
                        {
                            idxList.push_back(vSeq[i]);
                        }
                    }
                    int p = atoi(vItem->Attribute("ID"));
 
                    std::string vListStr =
                        ParseUtils::GenerateSeqString(idxList);
 
                    if (vListStr.length() == 0)
                    {
                        TiXmlElement *tmp = vItem;
                        vItem             = vItem->NextSiblingElement();
                        vBndRegions->RemoveChild(tmp);
                    }
                    else
                    {
                        vListStr                  = "C[" + vListStr + "]";
                        TiXmlText *vList          = new TiXmlText(vListStr);
                        TiXmlElement *vNewElement = new TiXmlElement("B");
                        vNewElement->SetAttribute("ID", p);
                        vNewElement->LinkEndChild(vList);
                        vNewBndRegions->LinkEndChild(vNewElement);
                        vBndRegionIdList.insert(p);
                        vItem = vItem->NextSiblingElement();
                    }
 
                    // store original bnd region order
                    m_bndRegOrder[p] = vSeq;
                }
                vConditions->ReplaceChild(vBndRegions, *vNewBndRegions);
            }
 
            if (vBndConditions)
            {
                vItem = vBndConditions->FirstChildElement();
                while (vItem)
                {
                    std::set<int>::iterator x;
                    if ((x = vBndRegionIdList.find(atoi(vItem->Attribute(
                             "REF")))) != vBndRegionIdList.end())
                    {
                        vItem->SetAttribute("REF", *x);
                        vItem = vItem->NextSiblingElement();
                    }
                    else
                    {
                        TiXmlElement *tmp = vItem;
                        vItem             = vItem->NextSiblingElement();
                        vBndConditions->RemoveChild(tmp);
                    }
                }
            }
            root->LinkEndChild(vConditions);
        }
 
        // Distribute other sections of the XML to each process as is.
        TiXmlElement *vSrc =
            m_session->GetElement("Nektar")->FirstChildElement();
        while (vSrc)
        {
            std::string vName = boost::to_upper_copy(vSrc->ValueStr());
            if (vName != "GEOMETRY" && vName != "CONDITIONS")
            {
                root->LinkEndChild(new TiXmlElement(*vSrc));
            }
            vSrc = vSrc->NextSiblingElement();
        }
 
        // Save Mesh
 
        boost::format pad("P%1$07d.xml");
        pad % partitions[i];
        boost::filesystem::path pFilename(pad.str());
 
        boost::filesystem::path fullpath = pdirname / pFilename;
        doc.SaveFile(LibUtilities::PortablePath(fullpath));
    }
}

References ASSERTL0, CellMLToNektar.pycml::format, CellMLToNektar.cellml_metadata::p, and Nektar::LibUtilities::PortablePath().

Member Data Documentation

className

std::string Nektar::SpatialDomains::MeshGraphXml::className

static

Initial value:

=
    GetMeshGraphFactory().RegisterCreatorFunction("Xml", MeshGraphXml::create,
                                                  "IO with Xml geometry")

Definition at line 67 of file MeshGraphXml.h.