gsFlowBndEvaluators.hpp

Go to the documentation of this file.

 
#pragma once
#include <gsIncompressibleFlow/src/gsFlowBndEvaluators.h>
 
namespace gismo
{
 
template<class T>
void gsFlowBndEvaluator<T>::evalOnPatchSide(index_t patchID, boxSide side)
{
    const gsBasis<T>* basis = &m_params.getBases().at(m_unkID).basis(patchID);
    short_t dim = m_params.getPde().domain().targetDim();
    m_mapData.patchId = patchID;
 
    gsVector<int> numQuadNodes(dim);
    const int dir = side.direction();
    for (short_t i = 0; i < dim; ++i)
        numQuadNodes[i] = basis->degree(i) + 1; // TODO: check correctness
    numQuadNodes[dir] = 1;
 
    gsGaussRule<T> QuRule(numQuadNodes);
 
    gsMatrix<T> quNodes; 
    gsVector<T> quWeights; 
 
    typename gsBasis<T>::domainIter domIt = basis->makeDomainIterator(side);
    for (; domIt->good(); domIt->next())
    {
        QuRule.mapTo(domIt->lowerCorner(), domIt->upperCorner(), quNodes, quWeights);
    
        m_mapData.points = quNodes;
        m_mapData.side = side;
        m_params.getPde().patches().patch(patchID).computeMap(m_mapData);
 
        this->evalOnElement(patchID, side, quNodes, quWeights);
    }
}
 
 
// ===================================================================================================================
 
 
template<class T>
void gsFlowBndEvaluator_flowRate<T>::evalOnElement(index_t patchID, boxSide side, const gsMatrix<T>& quNodes, const gsVector<T>& quWeights)
{
    gsMatrix<T> solUVals = m_velocityField.value(quNodes, patchID);
 
    for (index_t k = 0; k < quWeights.rows(); ++k)
    {
        gsVector<T> normal = m_mapData.outNormal(k);
        
        // the normal norm is equal to integral measure
        m_quantValue += quWeights[k] * normal.dot(solUVals.col(k));
    }
}
 
} // namespace gismo