This processing module calculates the Jacobian of elements using SpatialDomains::GeomFactors and the Element::GetGeom method. For now it simply prints a list of elements which have negative Jacobian. More...

#include <ProcessBL.h>

Inheritance diagram for Nektar::Utilities::ProcessBL:

[legend]

Collaboration diagram for Nektar::Utilities::ProcessBL:

[legend]

Public Member Functions
	ProcessBL (MeshSharedPtr m)
virtual	~ProcessBL ()
virtual void	Process ()
	Write mesh to output file.
Public Member Functions inherited from Nektar::Utilities::ProcessModule
	ProcessModule (FieldSharedPtr p_f)
	ProcessModule (MeshSharedPtr p_m)
Public Member Functions inherited from Nektar::Utilities::Module
	Module (FieldSharedPtr p_f)
virtual void	Process (po::variables_map &vm)=0
void	RegisterConfig (string key, string value)
	Register a configuration option with a module.
void	PrintConfig ()
	Print out all configuration options for a module.
void	SetDefaults ()
	Sets default configuration options for those which have not been set.
bool	GetRequireEquiSpaced (void)
void	SetRequireEquiSpaced (bool pVal)
	Module (MeshSharedPtr p_m)
void	RegisterConfig (string key, string value)
void	PrintConfig ()
void	SetDefaults ()
MeshSharedPtr	GetMesh ()
virtual void	ProcessVertices ()
	Extract element vertices.

Static Public Member Functions
static boost::shared_ptr< Module >	create (MeshSharedPtr m)
	Creates an instance of this class.

Static Public Attributes
static ModuleKey	className

Additional Inherited Members
Protected Member Functions inherited from Nektar::Utilities::Module
virtual void	ProcessEdges (bool ReprocessEdges=true)
	Extract element edges.
virtual void	ProcessFaces (bool ReprocessFaces=true)
	Extract element faces.
virtual void	ProcessElements ()
	Generate element IDs.
virtual void	ProcessComposites ()
	Generate composites.
void	ReorderPrisms (PerMap &perFaces)
	Reorder node IDs so that prisms and tetrahedra are aligned correctly.
void	PrismLines (int prism, PerMap &perFaces, set< int > &prismsDone, vector< ElementSharedPtr > &line)
Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr	m_f
	Field object.
map< string, ConfigOption >	m_config
	List of configuration values.
bool	m_requireEquiSpaced
MeshSharedPtr	m_mesh
	Mesh object.

Detailed Description

This processing module calculates the Jacobian of elements using SpatialDomains::GeomFactors and the Element::GetGeom method. For now it simply prints a list of elements which have negative Jacobian.

Definition at line 51 of file ProcessBL.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessBL::ProcessBL ( MeshSharedPtr m )

Definition at line 58 of file ProcessBL.cpp.

References Nektar::Utilities::Module::m_config.

                                            : ProcessModule(m)
        {
            // BL mesh configuration.
            m_config["layers"]     = ConfigOption(false, "2",
                "Number of layers to refine.");
            m_config["nq"]         = ConfigOption(false, "5",
                "Number of points in high order elements.");
            m_config["surf"]       = ConfigOption(false, "",
                "Tag identifying surface connected to prism.");
            m_config["r"]          = ConfigOption(false, "2.0",
                "Ratio to use in geometry progression.");
        }

Nektar::Utilities::ProcessBL::~ProcessBL ( )

virtual

Definition at line 71 of file ProcessBL.cpp.

        {
        }

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessBL::create ( MeshSharedPtr m )

inlinestatic

Creates an instance of this class.

Definition at line 55 of file ProcessBL.h.

                                                                   {
                return MemoryManager<ProcessBL>::AllocateSharedPtr(m);
            }

void Nektar::Utilities::ProcessBL::Process ( )

virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 76 of file ProcessBL.cpp.

        {
            if (m_mesh->m_verbose)
            {
                cout << "ProcessBL: Refining prismatic boundary layer..."
                     << endl;
            }
            int nodeId  = m_mesh->m_vertexSet.size();
            int nl      = m_config["layers"].as<int>();
            int nq      = m_config["nq"].    as<int>();
            // determine if geometric ratio is string or a constant. 
            LibUtilities::AnalyticExpressionEvaluator rEval;
            NekDouble r;
            int       rExprId       = -1;
            bool      ratioIsString = false;
            if (m_config["r"].isType<NekDouble>())
            {
                r = m_config["r"].as<NekDouble>();
            }
            else
            {
                std::string rstr = m_config["r"].as<string>();
                rExprId = rEval.DefineFunction("x y z", rstr);
                ratioIsString = true;
            }
            // Prismatic node -> face map.
            int prismFaceNodes[5][4] = {
                {0,1,2,3},{0,1,4,-1},{1,2,5,4},{3,2,5,-1},{0,3,5,4}};
            // Default PointsType.
            LibUtilities::PointsType pt = LibUtilities::eGaussLobattoLegendre;
            // Map which takes element ID to face on surface. This enables
            // splitting to occur in either y-direction of the prism.
            map<int, int> splitEls;
            // Set up maps which takes an edge (in nektar++ ordering) and return
            // their offset and stride in the 3d array of collapsed quadrature
            // points. Note that this map includes only the edges that are on
            // the triangular faces as the edges in the normal direction are
            // linear.
            int splitMapEdge[6] = {0,2,4,5,6,7};
            int splitMapOffset[6][2] = {
                {0,       1        },
                {nq,      1        },
                {0,       nq*(nl+1)},
                {nq-1,    nq*(nl+1)},
                {nq+nq-1, nq*(nl+1)},
                {nq,      nq*(nl+1)}
            };
            // splitEdge enumerates the edges in the standard prism along which
            // new nodes should be generated. These edges are the three between
            // the two triangular faces.
            int splitEdge[3] = {3,1,8};
            // edgeVertMap specifies the vertices which comprise those edges in
            // splitEdge; for example splitEdge[0] = 3 which connects vertices 0
            // and 3.
            int edgeVertMap[3][2] = {{0,3}, {1,2}, {4,5}};
            // edgeOffset holds the offset of each of edges 3, 1 and 8
            // respectively inside the collapsed coordinate system.
            int edgeOffset[3] = {0, nq-1, nq*(nl+1)*(nq-1)};
            // edgeMap associates geometry edge IDs to the (nl+1) vertices which
            // are generated along that edge when a prism is split, and is used
            // to avoid generation of duplicate vertices. It is stored as an
            // unordered map for speed.
            boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
            boost::unordered_map<int, vector<NodeSharedPtr> >::iterator eIt;
            string surf = m_config["surf"].as<string>();
            if (surf.size() > 0)
            {
                vector<unsigned int> surfs;
                ParseUtils::GenerateSeqVector(surf.c_str(), surfs);
                sort(surfs.begin(), surfs.end());
                // If surface is defined, process list of elements to find those
                // that are connected to it.
                for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
                {
                    ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
                    int nSurf = el->GetFaceCount();
                    for (int j = 0; j < nSurf; ++j)
                    {
                        int bl = el->GetBoundaryLink(j);
                        if (bl == -1)
                        {
                            continue;
                        }
                        ElementSharedPtr bEl  = m_mesh->m_element[m_mesh->m_expDim-1][bl];
                        vector<int>      tags = bEl->GetTagList();
                        vector<int>      inter;
                        sort(tags.begin(), tags.end());
                        set_intersection(surfs.begin(), surfs.end(),
                                         tags .begin(), tags .end(),
                                         back_inserter(inter));
                        ASSERTL0(inter.size() <= 1,
                                 "Intersection of surfaces wrong");
                        if (inter.size() == 1)
                        {
                            if (el->GetConf().m_e != LibUtilities::ePrism)
                            {
                                cerr << "WARNING: Found non-prismatic element "
                                     << "to split in surface " << surf
                                     << "; ignoring" << endl;
                                continue;
                            }
                            if (j % 2 == 0)
                            {
                                cerr << "WARNING: Found quadrilateral face "
                                     << j << " on surface " << surf
                                     << " connected to prism; ignoring."
                                     << endl;
                                continue;
                            }
                            if (splitEls.count(el->GetId()) > 0)
                            {
                                cerr << "WARNING: prism already found; "
                                     << "ignoring" << endl;
                            }
                            splitEls[el->GetId()] = j;
                        }
                    }
                }
            }
            else
            {
                // Otherwise, add all prismatic elements and assume face 1 of
                // the prism lies on the surface.
                for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
                {
                    ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
                    if (el->GetConf().m_e != LibUtilities::ePrism)
                    {
                        continue;
                    }
                    splitEls[el->GetId()] = 1;
                }
            }
            if (splitEls.size() == 0)
            {
                cerr << "WARNING: No elements detected to split." << endl;
                return;
            }
            // Erase all elements from the element list. Elements will be
            // re-added as they are split.
            vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
            m_mesh->m_element[m_mesh->m_expDim].clear();
            // Iterate over list of elements of expansion dimension.
            for (int i = 0; i < el.size(); ++i)
            {
                if (splitEls.count(el[i]->GetId()) == 0)
                {
                    m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
                    continue;
                }
                // Find quadrilateral boundary faces if any
                std::map<int, int> bLink;
                for (int j = 0; j < 5; j += 2)
                {
                    int bl = el[i]->GetBoundaryLink(j);
                    if (bl != -1)
                    {
                        bLink[j] = bl;
                    }
                }
                // Get elemental geometry object.
                SpatialDomains::PrismGeomSharedPtr geom =
                    boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(
                        el[i]->GetGeom(m_mesh->m_spaceDim));
                // Determine whether to use reverse points.
                LibUtilities::PointsType t = splitEls[el[i]->GetId()] == 1 ?
                    LibUtilities::eBoundaryLayerPoints :
                    LibUtilities::eBoundaryLayerPointsRev;
                if(ratioIsString) // determien value of r base on geom
                {
                    NekDouble x,y,z;
                    NekDouble x1,y1,z1;
                    int nverts = geom->GetNumVerts();
                    x = y = z = 0.0;
                    
                    for(int i = 0; i < nverts; ++i)
                    {
                        geom->GetVertex(i)->GetCoords(x1,y1,z1);
                        x += x1; y += y1; z += z1;
                    }
                    x /= (NekDouble) nverts;
                    y /= (NekDouble) nverts;
                    z /= (NekDouble) nverts;
                    r = rEval.Evaluate(rExprId,x,y,z,0.0);
                }
                // Create basis.
                LibUtilities::BasisKey B0(
                    LibUtilities::eModified_A, nq,
                    LibUtilities::PointsKey(nq,pt));
                LibUtilities::BasisKey B1(
                    LibUtilities::eModified_A, 2,
                    LibUtilities::PointsKey(nl+1, t, r));
                LibUtilities::BasisKey B2(
                    LibUtilities::eModified_B, nq,
                    LibUtilities::PointsKey(nq,pt));
                // Create local region.
                LocalRegions::PrismExpSharedPtr q =
                    MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
                        B0,B1,B2,geom);
                // Grab co-ordinates.
                Array<OneD, NekDouble> x(nq*nq*(nl+1));
                Array<OneD, NekDouble> y(nq*nq*(nl+1));
                Array<OneD, NekDouble> z(nq*nq*(nl+1));
                q->GetCoords(x,y,z);
                vector<vector<NodeSharedPtr> > edgeNodes(3);
                // Loop over edges to be split.
                for (int j = 0; j < 3; ++j)
                {
                    int locEdge = splitEdge[j];
                    int edgeId  = el[i]->GetEdge(locEdge)->m_id;
                    // Determine whether we have already generated vertices
                    // along this edge.
                    eIt = edgeMap.find(edgeId);
                    if (eIt == edgeMap.end())
                    {
                        // If not then resize storage to hold new points.
                        edgeNodes[j].resize(nl+1);
                        // Re-use existing vertices at endpoints of edge to
                        // avoid duplicating the existing vertices.
                        edgeNodes[j][0]  = el[i]->GetVertex(edgeVertMap[j][0]);
                        edgeNodes[j][nl] = el[i]->GetVertex(edgeVertMap[j][1]);
                        // Variable geometric ratio
                        if(ratioIsString) 
                        {
                            NekDouble x0,y0,z0;
                            NekDouble x1,y1,z1;
                            NekDouble xm,ym,zm;
                            
                            // -> Find edge end and mid points
                            x0 = x[edgeOffset[j]];
                            y0 = y[edgeOffset[j]];
                            z0 = z[edgeOffset[j]];
                            x1 = x[edgeOffset[j]+nl*nq];
                            y1 = y[edgeOffset[j]+nl*nq];
                            z1 = z[edgeOffset[j]+nl*nq];
                            xm = 0.5*(x0+x1);
                            ym = 0.5*(y0+y1);
                            zm = 0.5*(z0+z1);
                            // evaluate r factor based on mid point value
                            NekDouble rnew;
                            rnew = rEval.Evaluate(rExprId,xm,ym,zm,0.0);
                            // Get basis with new r; 
                            LibUtilities::PointsKey Pkey(nl+1, t, rnew);
                            LibUtilities::PointsSharedPtr newP
                                = LibUtilities::PointsManager()[Pkey];
                            
                            const Array<OneD, const NekDouble> z = newP->GetZ();
                            // Create new interior nodes based on this new blend
                            for (int k = 1; k < nl; ++k)
                            {
                                xm = 0.5*(1+z[k])*(x1-x0) + x0;
                                ym = 0.5*(1+z[k])*(y1-y0) + y0;
                                zm = 0.5*(1+z[k])*(z1-z0) + z0;
                                edgeNodes[j][k] = NodeSharedPtr(
                                        new Node(nodeId++, xm,ym,zm));
                            }
                        }
                        else
                        {
                            // Create new interior nodes.
                            for (int k = 1; k < nl; ++k)
                            {
                                int pos = edgeOffset[j] + k*nq;
                                edgeNodes[j][k] = NodeSharedPtr(
                                    new Node(nodeId++, x[pos], y[pos], z[pos]));
                            }
                        }
                        // Store these edges in edgeMap.
                        edgeMap[edgeId] = edgeNodes[j];
                    }
                    else
                    {
                        edgeNodes[j] = eIt->second;
                    }
                }
                // Create element layers.
                for (int j = 0; j < nl; ++j)
                {
                    // Offset of this layer within the collapsed coordinate
                    // system.
                    int offset = j*nq;
                    // Get corner vertices.
                    vector<NodeSharedPtr> nodeList(6);
                    nodeList[0] = edgeNodes[0][j  ];
                    nodeList[1] = edgeNodes[1][j  ];
                    nodeList[2] = edgeNodes[1][j+1];
                    nodeList[3] = edgeNodes[0][j+1];
                    nodeList[4] = edgeNodes[2][j  ];
                    nodeList[5] = edgeNodes[2][j+1];
                    // Create the element.
                    ElmtConfig conf(LibUtilities::ePrism, 1, true, true, false);
                    ElementSharedPtr elmt = GetElementFactory().
                        CreateInstance(
                      LibUtilities::ePrism,conf,nodeList,el[i]->GetTagList());
                    // Add high order nodes to split prismatic edges.
                    for (int l = 0; l < 6; ++l)
                    {
                        EdgeSharedPtr HOedge = elmt->GetEdge(
                            splitMapEdge[l]);
                        for (int k = 1; k < nq-1; ++k)
                        {
                            int pos = offset + splitMapOffset[l][0] +
                                k*splitMapOffset[l][1];
                            HOedge->m_edgeNodes.push_back(
                                NodeSharedPtr(
                                    new Node(nodeId++,x[pos],y[pos],z[pos])));
                        }
                        HOedge->m_curveType = pt;
                    }
                    // Change the surface elements of the quad face on the
                    // symmetry plane to match the layers of prisms.
                    map<int,int>::iterator it;
                    for (it = bLink.begin(); it != bLink.end(); ++it)
                    {
                        int fid = it->first;
                        int bl  = it->second;
                        if (j == 0)
                        {
                            // For first layer reuse existing 2D element.
                            ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim-1][bl];
                            for (int k = 0; k < 4; ++k)
                            {
                                e->SetVertex(
                                    k, nodeList[prismFaceNodes[fid][k]]);
                            }
                        }
                        else
                        {
                            // For all other layers create new element.
                            vector<NodeSharedPtr> qNodeList(4);
                            for (int k = 0; k < 4; ++k)
                            {
                                qNodeList[k] = nodeList[prismFaceNodes[fid][k]];
                            }
                            vector<int> tagBE;
                            tagBE = m_mesh->m_element[m_mesh->m_expDim-1][bl]->GetTagList();
                            ElmtConfig bconf(LibUtilities::eQuadrilateral,1,true,true,false);
                            ElementSharedPtr boundaryElmt = GetElementFactory().
                                CreateInstance(LibUtilities::eQuadrilateral,bconf,
                                               qNodeList,tagBE);
                            m_mesh->m_element[m_mesh->m_expDim-1].push_back(boundaryElmt);
                        }
                    }
                    m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
                }
            }
            // Re-process mesh to eliminate duplicate vertices and edges.
            ProcessVertices();
            ProcessEdges();
            ProcessFaces();
            ProcessElements();
            ProcessComposites();
        }

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessBL::className

static

Initial value:

            GetModuleFactory().RegisterCreatorFunction(
                ModuleKey(eProcessModule, "bl"), ProcessBL::create,
                "Refines a prismatic boundary layer.")

Definition at line 58 of file ProcessBL.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation