MultiphaseFlowComponent.cpp

#define _USE_MATH_DEFINES
 
#include <cmath>
#include <sstream>
#include "Types.h"
#include "Context.h"
#include "EntityStore.h"
 
#include "Utilities/Math.h"
#include "Systems/Core/CameraSystem.h"
#include "Systems/Core/MeshSystem.h"
#include "Math/CrossSectionGenerator.h"
 
#include "Resources/Mesh.h"
#include "Resources/MeshManager.h"
#include "Resources/MaterialManager.h"
#include "Resources/FontManager.h"
#include "Resources/DefaultMaterial.h"
 
#include "Components/Core/MeshRenderComponent.h"
#include "Components/Core/TextComponent.h"
#include "Components/Core/SceneComponent.h"
#include "Components/Core/MeshRenderComponent.h"
#include "Components/Core/TransformComponent.h"
#include "Components/Core/SubMeshRenderComponent.h"
#include "Components/Appearance/MaterialComponent.h"
 
#include "MultiphaseFlowComponent.h"
 
#include <glm/gtx/compatibility.hpp>
 
using namespace Cogs::Geometry;
using namespace Cogs::Reflection;
using namespace Cogs::Core;
 
static const int NUM_OF_SEGMENTS = 24;
static const int NUM_OF_SEGMENTS_SHADOW = 12;
static std::vector<glm::vec3> CROSS_SECTION;
static std::vector<glm::vec3> CROSS_SECTION_SHADOW;
 
static const glm::vec3 WATER_COLOR_DEFAULT = glm::vec3(54, 121, 159) / 255.f;
static const glm::vec3 OIL_COLOR_DEFAULT = glm::vec3(235, 201, 77) / 255.f;
static const glm::vec3 GAS_COLOR_DEFAULT = glm::vec3(195, 195, 195) / 255.f;
static const glm::vec3 SHADOW_COLOR = glm::vec3(0.3f, 0.3f, 0.3f);
static const glm::vec3 INDICATOR_RING_COLOR = glm::vec3(143, 14, 14) / 255.f;
 
static const float ANNULAR_FLOW_INCLINATION = 0.6f;
static const float GAS_STRATIFIED_TRANSPARENCY = 0.5f;
 
static const int GAS_STRATIFIED_TRANSITION_RANGE = 20;
 
MultiphaseFlowComponent::MultiphaseFlowComponent() :
    waterColor(WATER_COLOR_DEFAULT),
    oilColor(OIL_COLOR_DEFAULT),
    gasColor(GAS_COLOR_DEFAULT),
    indicatorRingColor(INDICATOR_RING_COLOR)
{
    if (CROSS_SECTION.empty())
    {
        CROSS_SECTION.resize(NUM_OF_SEGMENTS);
        CrossSectionGenerator::generateCircularCrossection(CROSS_SECTION.data(), static_cast<int>(NUM_OF_SEGMENTS), 1.f, 0, glm::pi<float>() * 2.0f, false);
    }
 
    if (CROSS_SECTION_SHADOW.empty())
    {
        CROSS_SECTION_SHADOW.resize(NUM_OF_SEGMENTS_SHADOW);
        CrossSectionGenerator::generateCircularCrossection(CROSS_SECTION_SHADOW.data(), static_cast<int>(NUM_OF_SEGMENTS_SHADOW), 1.f, 0, glm::pi<float>() * 2.0f, false);
    }
}
 
void MultiphaseFlowComponent::registerType()
{
    Field fields[] = {
        Field(Name("trajectoryPoints"), &MultiphaseFlowComponent::trajectoryPoints),
        Field(Name("radius"), &MultiphaseFlowComponent::radius),
        Field(Name("flowFraction"), &MultiphaseFlowComponent::flowFraction),
        Field(Name("wallValues"), &MultiphaseFlowComponent::wallValues),
        Field(Name("heightScale"), &MultiphaseFlowComponent::heightScale),
        Field(Name("radiusScale"), &MultiphaseFlowComponent::radiusScale),
        Field(Name("showWater"), &MultiphaseFlowComponent::showWater),
        Field(Name("showOil"), &MultiphaseFlowComponent::showOil),
        Field(Name("showGas"), &MultiphaseFlowComponent::showGas),
        Field(Name("unitScale"), &MultiphaseFlowComponent::unitScale),
        Field(Name("unitText"), &MultiphaseFlowComponent::unitText),
        Field(Name("updateNeeded"), &MultiphaseFlowComponent::updateNeeded),
        Field(Name("waterColor"), &MultiphaseFlowComponent::waterColor),
        Field(Name("oilColor"), &MultiphaseFlowComponent::oilColor),
        Field(Name("gasColor"), &MultiphaseFlowComponent::gasColor),
        Field(Name("fontSize"), &MultiphaseFlowComponent::fontSize),
        Field(Name("showAxialExtents"), &MultiphaseFlowComponent::showAxialExtents),
        Field(Name("showShadow"), &MultiphaseFlowComponent::showShadow),
        Field(Name("showGrid"), &MultiphaseFlowComponent::showGrid),
        Field(Name("centerMesh"), &MultiphaseFlowComponent::centerMesh),
        Field(Name("showIndicatorRing"), &MultiphaseFlowComponent::showIndicatorRing),
        Field(Name("indicatorRingPosition"), &MultiphaseFlowComponent::indicatorRingPosition)
    };
 
    Method methods[] = {
        Method(Name("initialize"), &MultiphaseFlowComponent::initialize),
        Method(Name("update"), &MultiphaseFlowComponent::update),
    };
 
    DynamicComponent::registerDerivedType<MultiphaseFlowComponent>().setFields(fields).setMethods(methods);
}
 
bool MultiphaseFlowComponent::isStratified(MeshType liquidType)
{
    assert(liquidType != MeshType::Wall && liquidType != MeshType::Shadow);
 
    return (liquidType == MeshType::StratifiedWater ||
        liquidType == MeshType::StratifiedOil ||
        liquidType == MeshType::StratifiedGas);
}
 
bool MultiphaseFlowComponent::isAnnular(MeshType liquidType)
{
    assert(liquidType != MeshType::Wall && liquidType != MeshType::Shadow);
 
    return !isStratified(liquidType);
}
 
unsigned int MultiphaseFlowComponent::getNumOfSegmentsInCrossSection(MeshType liquidType)
{
    return liquidType != MeshType::Shadow ? NUM_OF_SEGMENTS : NUM_OF_SEGMENTS_SHADOW;
}
 
void MultiphaseFlowComponent::createMeshEntity(Cogs::Core::EntityPtr& mesh, MeshType liquidType)
{
    mesh = context->store->createChildEntity("MeshPart", getContainer());
 
    MeshComponent* meshComp = mesh->getComponent<MeshComponent>();
    if (meshComp->meshHandle == MeshHandle::NoHandle)
    {
        meshComp->meshHandle = context->meshManager->create();
    }
 
    auto materialInstance = context->materialInstanceManager->createMaterialInstance(MaterialHandle(mainMaterial));
    auto meshRenderComp = mesh->getComponent<MeshRenderComponent>();
    meshRenderComp->material = materialInstance;
 
    updateMaterial(mesh, liquidType);
}
 
void MultiphaseFlowComponent::updateMaterial(Cogs::Core::EntityPtr mesh, MeshType liquidType)
{
    // Set default cull mode.
    auto meshRenderComp = mesh->getComponent<MeshRenderComponent>();
    auto material = meshRenderComp->material;
    material->setOption("CullMode", "Front");
 
    // Set material color and transparency based on the liquid type.
    glm::vec3 diffuseColor;
    switch (liquidType)
    {
    case MeshType::StratifiedWater:
        diffuseColor = waterColor;
        material->setOption("CullMode", "Back");
        material->setOption("DepthBiasEnabled", "true");
        material->setOption("DepthBiasConstant", "0.5f");
        material->setOption("DepthBiasSlope", "0.5f");
        break;
    case MeshType::AnnularWater:
        diffuseColor = waterColor;
        material->setOption("CullMode", "Back");
        break;
    case MeshType::StratifiedOil:
        diffuseColor = oilColor;
        material->setOption("CullMode", "None");
        break;
    case MeshType::AnnularOil:
        diffuseColor = oilColor;
        material->setOption("Transparency", "On");
        break;
    case MeshType::StratifiedGas:
        diffuseColor = gasColor;
        material->setOption("Transparency", "On");
        break;
    case MeshType::AnnularGas:
        diffuseColor = gasColor;
        break;
    case MeshType::Wall:
        diffuseColor = glm::vec3(1.f);
        material->setVariant("EnablePerVertexColor", true);
        break;
    case MeshType::Shadow:
        diffuseColor = SHADOW_COLOR;
        material->setVariant("EnableLighting", false);
        material->setVariant("LightModel", "BaseColor");
        material->setOption("CullMode", "None");
        break;
    default:
        assert(false && "Unknown liquid type");
        break;
    }
 
    // Set material ambient and specular color.
    material->setVec3Property(diffuseColorMaterialPropertyKey, diffuseColor);
    material->setVec3Property(emissiveColorMaterialPropertyKey, diffuseColor * (liquidType != MeshType::Wall ? 0.3f : 0.05f));
    material->setVec3Property(specularColorMaterialPropertyKey, liquidType != MeshType::Shadow ? diffuseColor : glm::vec3(0, 0, 0));
}
 
void MultiphaseFlowComponent::initialize(Context * ctx)
{
    context = ctx;
 
    // Create material.
    auto materialHandle = context->materialManager->loadMaterial("MultiphaseFlowMaterial.material");
    materialHandle->setMaterialFlag(MaterialFlags::Default);
    context->materialManager->processLoading();
    materialHandle->setChanged();
    mainMaterial = materialHandle.resolve();
 
    diffuseColorMaterialPropertyKey = materialHandle->getVec3Key("diffuseColor");
    emissiveColorMaterialPropertyKey = materialHandle->getVec3Key("emissiveColor");
    specularColorMaterialPropertyKey = materialHandle->getVec3Key("specularColor");
 
    // Initialize meshes for all liquids.
    for (int i = MeshType::First; i <= MeshType::Last; i++)
    {
        auto liquidType = MeshType(i);
        createMeshEntity(meshEntity[i], liquidType);
    }
 
    // Create ground mesh entity.
    createGridEntity();
 
    // Create indicator ring.
    createIndicatorRing();
 
    // Create text entity.
    createTextEntity();
}
 
static void SetLineMaterial(Context* context, EntityPtr entity, const glm::vec3& color, float lineWidth = 1.f)
{
    auto material = context->materialInstanceManager->createMaterialInstance(context->materialManager->getDefaultMaterial());
    material->setPermutation("Line");
    material->setVariant("EnableLighting", false);
    material->setBoolProperty(DefaultMaterial::EnableLighting, false);
    material->setFloatProperty(DefaultMaterial::LineWidth, lineWidth);
    material->setVec4Property(DefaultMaterial::DiffuseColor, glm::vec4(color, 1));
 
    auto meshRenderComp = entity->getComponent<MeshRenderComponent>();
    meshRenderComp->material = material;
}
 
// Create ground mesh entity.
void MultiphaseFlowComponent::createGridEntity()
{
    groundGridMeshEntity = context->store->createChildEntity("MeshPart", getContainer());
    auto* meshComp = groundGridMeshEntity->getComponent<MeshComponent>();
    if (meshComp->meshHandle == MeshHandle::NoHandle)
    {
        meshComp->meshHandle = context->meshManager->create();
    }
    SetLineMaterial(context, groundGridMeshEntity, glm::vec3(0.6f, 0.6f, 0.6f));
 
    annotationAxisMeshEntity = context->store->createChildEntity("MeshPart", getContainer());
    auto* annotationMeshComp = annotationAxisMeshEntity->getComponent<MeshComponent>();
    if (annotationMeshComp->meshHandle == MeshHandle::NoHandle)
    {
        annotationMeshComp->meshHandle = context->meshManager->create();
    }
    SetLineMaterial(context, annotationAxisMeshEntity, glm::vec3(0.2f, 0.2f, 0.2f));
}
 
void MultiphaseFlowComponent::createIndicatorRing()
{
    indicatorRingEntity = context->store->createChildEntity("MeshPart", getContainer());
    auto* meshComp = indicatorRingEntity->getComponent<MeshComponent>();
    if (meshComp->meshHandle == MeshHandle::NoHandle)
    {
        meshComp->meshHandle = context->meshManager->create();
    }
 
    SetLineMaterial(context, indicatorRingEntity, indicatorRingColor, 5.f);
}
 
// Create text entity.
void MultiphaseFlowComponent::createTextEntity()
{
    for (int i = AnnotationAxis::AxisFirst; i <= AnnotationAxis::AxisLast; i++)
    {
        textEntity[i] = context->store->createChildEntity("Text", getContainer());
        auto* textComp = textEntity[i]->getComponent<TextComponent>();
        textComp->positionMode = PositionMode::World;
        textComp->horizontalJustification = HorizontalJustification::Left;
        textComp->verticalAlignment = i == AnnotationAxis::X ? VerticalAlignment::Top : VerticalAlignment::Center;
    }
 
    assert(currentFontSize != fontSize);
    updateFontSize();
}
 
void MultiphaseFlowComponent::updateFontSize()
{
    if (currentFontSize != fontSize)
    {
        currentFontSize = fontSize;
 
        for (int i = AnnotationAxis::AxisFirst; i <= AnnotationAxis::AxisLast; i++)
        {
            auto* textComp = textEntity[i]->getComponent<TextComponent>();
            textComp->fontHandle = context->fontManager->loadFont("Fonts/SourceSansPro-Regular.ttf", currentFontSize, NoResourceId);
        }
    }
}
 
static float sCurve(float a_fValue, float a_fPower = 1.5f)
{
    if (a_fValue < 0.5f)
    {
        return (glm::pow(glm::clamp(a_fValue, 0.f, 1.f) * 2, a_fPower)) * 0.5f;
    }
    else
    {
        return -glm::pow(-(glm::clamp(a_fValue, 0.f, 1.f) - 1) * 2, a_fPower) * 0.5f + 1;
    }
}
 
static glm::vec3 wavelengthToColor(float lambda)
{
    // Based on: http://www.efg2.com/Lab/ScienceAndEngineering/Spectra.htm
    // The foregoing is based on: http://www.midnightkite.com/color.html
    struct Color
    {
        float c[3];
        glm::vec3 toColor(float factor) const
        {
            float const gamma = 0.8f;
            float ci[3];
            for (int i = 0; i < 3; ++i)
            {
                ci[i] = c[i] == 0 ? 0 : glm::round(255 * pow(c[i] * factor, gamma));
            }
            return glm::vec3(ci[0] / 255.f, ci[1] / 255.f, ci[2] / 255.f);
        }
    } color;
 
    float factor = 0.f;
 
    color.c[0] = 0.f;
    color.c[1] = 0.f;
    color.c[2] = 0.f;
 
    static float thresholds[] = { 380, 440, 490, 510, 580, 645, 780 };
 
    for (unsigned int i = 0; i < sizeof(thresholds) / sizeof(thresholds[0]); ++i)
    {
        float t1 = thresholds[i], t2 = thresholds[i + 1];
        if (lambda < t1 || lambda >= t2)
        {
            continue;
        }
        if (i % 2)
        {
            std::swap(t1, t2);
        }
        color.c[i % 3] = (i < 5) ? (lambda - t2) / (t1 - t2) : 0.f;
        color.c[2 - i / 2] = 1.f;
        factor = 1.f;
        break;
    }
 
    // Let the intensity fall off near the vision limits
    if (lambda >= 380 && lambda < 420)
    {
        factor = 0.3f + 0.7f * (lambda - 380) / (420 - 380);
    }
    else if (lambda >= 700 && lambda < 780)
    {
        factor = 0.3f + 0.7f * (780 - lambda) / (780 - 700);
    }
    return color.toColor(factor);
}
 
glm::vec3 MultiphaseFlowComponent::wallValueToColor(const ControlPoint& controlPoint)
{
    static const float f1 = 1.f / 435.f;
    static const float f2 = 1.f / 650.f;
    float lambda = 1.f / (f1 - controlPoint.wall * (f1 - f2));
    return wavelengthToColor(lambda);
}
 
void MultiphaseFlowComponent::transformCrossection(float annularRatio, float alpha, std::vector<Vertex>& vertexBuffer, const ControlPoint& controlPoint, MeshType liquidType, bool isClosingCap /*= false*/)
{
    // Modify radius.
    float radius = controlPoint.radius;
    float annularOilTransparency = 1.f;
    if (liquidType == MeshType::AnnularOil && !isClosingCap)
    {
        float radiusAnnular = glm::sqrt(controlPoint.flow.oil + controlPoint.flow.gas);
        radius *= glm::lerp(1.0f, radiusAnnular, annularRatio);
 
        annularOilTransparency = glm::lerp(1.f, controlPoint.flow.oil, annularRatio);
    }
    else
        if (liquidType == MeshType::AnnularGas)
        {
            float radiusAnnular = glm::sqrt(controlPoint.flow.gas);
            radius *= glm::lerp(1.0f, radiusAnnular, annularRatio);
        }
 
    // Reduce multiphase flow mesh radius when the wall is visible.
    bool wallVisible = (wallValues.size() == trajectoryPoints.size());
    if (wallVisible && liquidType != MeshType::Wall && liquidType != MeshType::Shadow)
    {
        radius *= 0.8f;
    }
 
    // Process all crosssection points
    const auto& crossSection = (liquidType != MeshType::Shadow ? CROSS_SECTION : CROSS_SECTION_SHADOW);
    for (size_t j = 0; j < crossSection.size(); ++j)
    {
        glm::vec3 positionN = controlPoint.rotation * crossSection[j];
        glm::vec3 crossectionPointAnnular = crossSection[j];
 
        // Modify position based on the liquid fraction for stratified flow.
        glm::vec3 crossectionPointStratified = crossSection[j];
        glm::vec3 crossectionPointStratifiedInverted = crossSection[j];
        switch (liquidType)
        {
        case MeshType::AnnularWater:
            break;
 
        case MeshType::StratifiedWater:
 
            if (crossectionPointStratified.x < controlPoint.flow.heightWaterOil)
            {
                crossectionPointStratified.x = controlPoint.flow.heightWaterOil;
                crossectionPointStratified.y = crossectionPointStratified.y > 0.f ? controlPoint.flow.widthWaterOil : -controlPoint.flow.widthWaterOil;
            }
 
            // Inverted.
            if (crossectionPointStratifiedInverted.x > -controlPoint.flow.heightWaterOil)
            {
                crossectionPointStratifiedInverted.x = -controlPoint.flow.heightWaterOil;
                crossectionPointStratifiedInverted.y = crossectionPointStratifiedInverted.y > 0.f ? controlPoint.flow.widthWaterOil : -controlPoint.flow.widthWaterOil;
            }
            break;
 
        case MeshType::StratifiedOil:
        case MeshType::AnnularOil:
            if (crossectionPointStratified.x < controlPoint.flow.heightOilGas)
            {
                crossectionPointStratified.x = controlPoint.flow.heightOilGas;
                crossectionPointStratified.y = crossectionPointStratified.y > 0.f ? controlPoint.flow.widthOilGas : -controlPoint.flow.widthOilGas;
            }
            else
                if (crossectionPointStratified.x > controlPoint.flow.heightWaterOil)
                {
                    crossectionPointStratified.x = controlPoint.flow.heightWaterOil;
                    crossectionPointStratified.y = crossectionPointStratified.y > 0.f ? controlPoint.flow.widthWaterOil : -controlPoint.flow.widthWaterOil;
                }
 
            // Inverted.
            if (crossectionPointStratifiedInverted.x > -controlPoint.flow.heightOilGas)
            {
                crossectionPointStratifiedInverted.x = -controlPoint.flow.heightOilGas;
                crossectionPointStratifiedInverted.y = crossectionPointStratifiedInverted.y > 0.f ? controlPoint.flow.widthOilGas : -controlPoint.flow.widthOilGas;
            }
            else
                if (crossectionPointStratifiedInverted.x < -controlPoint.flow.heightWaterOil)
                {
                    crossectionPointStratifiedInverted.x = -controlPoint.flow.heightWaterOil;
                    crossectionPointStratifiedInverted.y = crossectionPointStratifiedInverted.y > 0.f ? controlPoint.flow.widthWaterOil : -controlPoint.flow.widthWaterOil;
                }
            break;
 
        case MeshType::StratifiedGas:
        case MeshType::AnnularGas:
            if (crossectionPointStratified.x > controlPoint.flow.heightOilGas)
            {
                crossectionPointStratified.x = controlPoint.flow.heightOilGas;
                crossectionPointStratified.y = crossectionPointStratified.y > 0.f ? controlPoint.flow.widthOilGas : -controlPoint.flow.widthOilGas;
            }
 
            // Inverted.
            if (crossectionPointStratifiedInverted.x < -controlPoint.flow.heightOilGas)
            {
                crossectionPointStratifiedInverted.x = -controlPoint.flow.heightOilGas;
                crossectionPointStratifiedInverted.y = crossectionPointStratifiedInverted.y > 0.f ? controlPoint.flow.widthOilGas : -controlPoint.flow.widthOilGas;
            }
            break;
 
        case MeshType::Shadow:
        case MeshType::Wall:
            break;
 
        default:
            break;
        }
 
        // Interpolate between inverted and non-inverted cross sections (oil always on the top).
        if (liquidType == MeshType::StratifiedWater)
        {
            // Special treatment for stratified oil - binary interpolation.
            crossectionPointStratified = (controlPoint.invertedFlowRatio < 0.5f ? crossectionPointStratified : crossectionPointStratifiedInverted);
        }
        else
            if (liquidType != MeshType::Shadow && liquidType != MeshType::Wall)
            {
                crossectionPointStratified = glm::lerp(crossectionPointStratified, crossectionPointStratifiedInverted, controlPoint.invertedFlowRatio);
            }
 
        // Interpolate between stratified and annular flow.
        glm::vec3 crossectionPoint = glm::lerp(crossectionPointStratified, crossectionPointAnnular, annularRatio);
 
        // Calculate rotated position and normal vector.
        glm::vec3 position = controlPoint.rotation * crossectionPoint;
 
        // Set radius and apply offset.
        position *= radius;
        position += controlPoint.position;
 
        // Create vertex color.
        auto transparency = 1.f;
        switch (liquidType)
        {
        case MeshType::StratifiedWater:
            break;
        case MeshType::AnnularWater:
            break;
        case MeshType::StratifiedOil:
            break;
        case MeshType::AnnularOil:
            transparency = annularOilTransparency;
            break;
        case MeshType::StratifiedGas:
            transparency = alpha;
            break;
        case MeshType::AnnularGas:
            break;
        case MeshType::Wall:
            break;
        case MeshType::Shadow:
            break;
        default:
            assert(false && "Unknown liquid type");
            break;
        }
 
        // Get normal.
        glm::vec3 normal;
        if (!isClosingCap)
        {
            normal = glm::normalize(positionN);
        }
        else
        {
            // Normal vector for closing caps based on the the trajectory direction.
            normal = controlPoint.direction;
        }
 
        // Construct final vertex.
        vertexBuffer.push_back({ position, glm::vec4(-normal, transparency) });
    }
}
 
void MultiphaseFlowComponent::generateExtrusionVertices(std::vector<Vertex>& vertexBuffer, std::vector<unsigned int>& indexBuffer, const std::vector<ControlPoint>& controlPoints, MeshType liquidType)
{
    assert(controlPoints.size() >= 2);
 
    bool isAnnular = MultiphaseFlowComponent::isAnnular(liquidType);
 
    // Reserve memory for the vertex buffer.
    assert(vertexBuffer.size() == 0);
    vertexBuffer.reserve(NUM_OF_SEGMENTS * controlPoints.size() + NUM_OF_SEGMENTS * 4); // Allocate additional memory for closing caps.
 
                                                                                        // Create flow ranges.
    std::vector<std::pair<int, int>> flowRanges;
    int startRange = -1;
    for (size_t i = 0; i < controlPoints.size(); i++)
    {
        bool isPointAnnularFlow = (controlPoints[i].annularFlowRatio > 0.f);
        if (isPointAnnularFlow == isAnnular)
        {
            if (startRange == -1)
            {
                startRange = (int)i;
            }
        }
        else
        {
            if (startRange >= 0)
            {
                assert(i > 0);
                flowRanges.push_back(std::pair<int, int>(startRange, (int)i));
                startRange = -1;
            }
        }
    }
 
    if (startRange >= 0)
    {
        flowRanges.push_back(std::pair<int, int>(startRange, (int)controlPoints.size() - 1));
    }
 
    // Generate geometry for all ranges.
    unsigned int segmentStart = 0;
    for (const auto& range : flowRanges)
    {
        int segmentRange = range.second - range.first;
        if (segmentRange <= 0)
        {
            continue;
        }
 
        // Fill up vertex buffer with transformed vertices.
        for (int i = range.first; i <= range.second; i++)
        {
            const auto& controlPoint = controlPoints[i];
            float annularRatio = controlPoint.annularFlowRatio;
            float alpha = 1.f;
 
            if ((i == range.first || i == range.second) && i > 0 && i < static_cast<int>(controlPoints.size()) - 1)
            {
                annularRatio = 0.f;
            }
 
            // Special transition for stratified gas.
            if (liquidType == MeshType::StratifiedGas)
            {
                float ratioNext = 0;
                for (int j = 0; j <= GAS_STRATIFIED_TRANSITION_RANGE; j++)
                {
                    int nextIndex = i + j;
                    if (nextIndex >= static_cast<int>(controlPoints.size()))
                    {
                        break;
                    }
                    const auto& nextCtrlPoint = controlPoints[nextIndex];
 
                    if (nextCtrlPoint.annularFlowRatio > 0.f)
                    {
                        ratioNext = 1.f - static_cast<float>(j) / static_cast<float>(GAS_STRATIFIED_TRANSITION_RANGE);
                        break;
                    }
                }
 
                float ratioPrev = 0;
                for (int j = 0; j <= GAS_STRATIFIED_TRANSITION_RANGE; j++)
                {
                    int prevIndex = i - j;
                    if (prevIndex < 0)
                    {
                        break;
                    }
                    const auto& nextCtrlPoint = controlPoints[prevIndex];
 
                    if (nextCtrlPoint.annularFlowRatio > 0.f)
                    {
                        ratioPrev = 1.f - static_cast<float>(j) / static_cast<float>(GAS_STRATIFIED_TRANSITION_RANGE);
                        break;
                    }
                }
 
                alpha = std::max(ratioNext, ratioPrev);
                alpha = GAS_STRATIFIED_TRANSPARENCY + alpha * (1.f - GAS_STRATIFIED_TRANSPARENCY);
            }
 
            transformCrossection(annularRatio, alpha, vertexBuffer, controlPoint, liquidType);
        }
 
        assert(segmentRange > 0);
        unsigned int segmentEnd = segmentStart + segmentRange;
        assert(segmentStart < segmentEnd && segmentEnd);
 
        // Generate index buffer.
        for (unsigned int offset = segmentStart; offset <= segmentEnd - 1; offset++)
        {
            const auto numOfSegments = getNumOfSegmentsInCrossSection(liquidType);
            for (unsigned int j = 0; j < numOfSegments - 1; j++)
            {
                indexBuffer.push_back((offset * numOfSegments) + j);
                indexBuffer.push_back((offset * numOfSegments) + j + 1);
                indexBuffer.push_back(((offset + 1) * numOfSegments) + j + 1);
 
                indexBuffer.push_back((offset * numOfSegments) + j);
                indexBuffer.push_back(((offset + 1) * numOfSegments) + j + 1);
                indexBuffer.push_back(((offset + 1) * numOfSegments) + j);
            }
 
            indexBuffer.push_back((offset * numOfSegments) + numOfSegments - 1);
            indexBuffer.push_back((offset * numOfSegments));
            indexBuffer.push_back(((offset + 1) * numOfSegments));
 
            indexBuffer.push_back((offset * numOfSegments) + numOfSegments - 1);
            indexBuffer.push_back(((offset + 1) * numOfSegments));
            indexBuffer.push_back(((offset + 1) * numOfSegments) + numOfSegments - 1);
        }
 
        segmentStart = segmentEnd + 1;
    }
 
    // Closing caps for solid meshes (stratified oil and annular gas).
    if (liquidType == MeshType::StratifiedOil ||
        liquidType == MeshType::AnnularGas)
    {
        // For annular gas flow crate closing caps at bot ends of every range.
        for (const auto& range : flowRanges)
        {
            int segmentRange = range.second - range.first;
            if (segmentRange <= 0)
            {
                continue;
            }
 
            if (range.first > 0)
            {
                generateClosingCaps(vertexBuffer, indexBuffer, controlPoints[range.first], liquidType);
            }
 
            if (range.second < static_cast<int>(controlPoints.size()) - 1)
            {
                generateClosingCaps(vertexBuffer, indexBuffer, controlPoints[range.second], liquidType);
            }
        }
    }
 
    // Create closing caps at the both ends of the flow.
    if (!flowRanges.empty() && liquidType != MeshType::Wall && liquidType != MeshType::Shadow)
    {
        if (flowRanges[0].first == 0)
        {
            generateClosingCaps(vertexBuffer, indexBuffer, controlPoints[0], liquidType);
        }
 
        if (flowRanges[flowRanges.size() - 1].second == static_cast<int>(controlPoints.size()) - 1)
        {
            generateClosingCaps(vertexBuffer, indexBuffer, controlPoints[controlPoints.size() - 1], liquidType);
        }
    }
}
 
void MultiphaseFlowComponent::generateClosingCaps(std::vector<Vertex>& vertexBuffer, std::vector<unsigned int>& indexBuffer, const ControlPoint& controlPoint, MeshType liquidType)
{
    float ratio = (controlPoint.annularFlowRatio < 0.5f ? 0.f : 1.f);
    float alpha = 1.f;
 
    if (liquidType == MeshType::StratifiedGas)
    {
        alpha = GAS_STRATIFIED_TRANSPARENCY;
    }
 
    if (liquidType == MeshType::Wall || isStratified(liquidType) || liquidType == MeshType::AnnularGas)
    {
        for (int i = 0; i < 2; i++)
        {
            ControlPoint ctrlPoint = controlPoint;
 
            if (i == 0)
            {
                ctrlPoint.direction *= -1;
            }
 
            const auto vertexBufferOffset = static_cast<unsigned int>(vertexBuffer.size());
 
            // Generate vertex buffer for closing caps.
            transformCrossection(ratio, alpha, vertexBuffer, ctrlPoint, liquidType, true);
 
            const auto numOfSegments = getNumOfSegmentsInCrossSection(liquidType);
            for (unsigned int j = 0; j < numOfSegments - 2; ++j)
            {
                indexBuffer.push_back(vertexBufferOffset);
                if (i == 1)
                {
                    indexBuffer.push_back(vertexBufferOffset + j + 2);
                    indexBuffer.push_back(vertexBufferOffset + j + 1);
                }
                else
                {
                    indexBuffer.push_back(vertexBufferOffset + j + 1);
                    indexBuffer.push_back(vertexBufferOffset + j + 2);
                }
            }
        }
    }
    else
        if (liquidType == MeshType::AnnularWater || liquidType == MeshType::AnnularOil)
        {
            assert(isAnnular(liquidType));
 
            const unsigned int offset = static_cast<unsigned int>(vertexBuffer.size());
 
            ControlPoint ctrlPoint1 = controlPoint;
            ControlPoint ctrlPoint2 = controlPoint;
 
            if (liquidType == MeshType::AnnularWater)
            {
                ctrlPoint2.radius *= glm::sqrt(ctrlPoint2.flow.oil + ctrlPoint2.flow.gas);
 
                ctrlPoint1.direction *= -1;
                ctrlPoint2.direction *= -1;
            }
            else
            {
                ctrlPoint1.radius *= glm::sqrt(ctrlPoint1.flow.oil + ctrlPoint1.flow.gas);
                ctrlPoint2.radius *= glm::sqrt(ctrlPoint2.flow.gas);
            }
 
            // Generate vertex buffer for closing caps.
            transformCrossection(ratio, alpha, vertexBuffer, ctrlPoint1, liquidType, true);
            transformCrossection(ratio, alpha, vertexBuffer, ctrlPoint2, liquidType, true);
 
            const auto numOfSegments = getNumOfSegmentsInCrossSection(liquidType);
            for (unsigned int i = 0; i < 2; i++)
            {
                for (unsigned int j = 0; j < numOfSegments; ++j)
                {
                    indexBuffer.push_back(offset + j);
                    if (i == 0)
                    {
                        indexBuffer.push_back(offset + j + numOfSegments);
                        indexBuffer.push_back(j < numOfSegments - 1 ? offset + j + 1 : offset);
                    }
                    else
                    {
                        indexBuffer.push_back(j < numOfSegments - 1 ? offset + j + 1 : offset);
                        indexBuffer.push_back(offset + j + numOfSegments);
                    }
 
                    indexBuffer.push_back(offset + j + numOfSegments);
                    if (i == 0)
                    {
                        indexBuffer.push_back(j < numOfSegments - 1 ? offset + j + numOfSegments + 1 : offset + numOfSegments);
                        indexBuffer.push_back(j < numOfSegments - 1 ? offset + j + 1 : offset);
                    }
                    else
                    {
                        indexBuffer.push_back(j < numOfSegments - 1 ? offset + j + 1 : offset);
                        indexBuffer.push_back(j < numOfSegments - 1 ? offset + j + numOfSegments + 1 : offset + numOfSegments);
                    }
                }
            }
        }
        else
        {
            assert(!"Unknown flow type");
        }
}
 
void MultiphaseFlowComponent::generateExtrusionVerticesWall(std::vector<WallVertex>& vertexBuffer, std::vector<unsigned int>& indexBuffer, const std::vector<ControlPoint>& controlPoints, MeshType liquidType)
{
    for (const auto& controlPoint : controlPoints)
    {
        std::vector<Vertex> temp;
        transformCrossection(1.f, 1.f, temp, controlPoint, liquidType);
 
        // Add colors.
        for (const auto& v : temp)
        {
            vertexBuffer.push_back({ v.position, v.normalAndOpacity, wallValueToColor(controlPoint) });
        }
    }
 
    // Generate index buffer.
    const auto numOfSegments = getNumOfSegmentsInCrossSection(liquidType);
    for (unsigned int offset = 0; offset < static_cast<unsigned int>(controlPoints.size()) - 1; offset++)
    {
        for (unsigned int j = 0; j < numOfSegments - 1; j++)
        {
            indexBuffer.push_back((offset * numOfSegments) + j);
            indexBuffer.push_back((offset * numOfSegments) + j + 1);
            indexBuffer.push_back(((offset + 1) * numOfSegments) + j + 1);
 
            indexBuffer.push_back((offset * numOfSegments) + j);
            indexBuffer.push_back(((offset + 1) * numOfSegments) + j + 1);
            indexBuffer.push_back(((offset + 1) * numOfSegments) + j);
        }
 
        indexBuffer.push_back((offset * numOfSegments) + numOfSegments - 1);
        indexBuffer.push_back((offset * numOfSegments));
        indexBuffer.push_back(((offset + 1) * numOfSegments));
 
        indexBuffer.push_back((offset * numOfSegments) + numOfSegments - 1);
        indexBuffer.push_back(((offset + 1) * numOfSegments));
        indexBuffer.push_back(((offset + 1) * numOfSegments) + numOfSegments - 1);
    }
 
    // Generate closing caps if multiphase flow mesh is not visible.
    if (liquidType == MeshType::Wall && flowFraction.size() != trajectoryPoints.size())
    {
        std::vector<Vertex> tempFirst;
        std::vector<Vertex> tempLast;
        std::vector<unsigned int> indexBufferFirst;
        std::vector<unsigned int> indexBufferLast;
        generateClosingCaps(tempFirst, indexBufferFirst, controlPoints[0], liquidType);
        generateClosingCaps(tempLast, indexBufferLast, controlPoints[controlPoints.size() - 1], liquidType);
 
        // Add colors.
        unsigned int offsetFirst = static_cast<unsigned int>(vertexBuffer.size());
        vertexBuffer.reserve(vertexBuffer.size() + tempFirst.size() + tempLast.size());
        for (const auto& v : tempFirst)
        {
            vertexBuffer.push_back({ v.position, v.normalAndOpacity, wallValueToColor(controlPoints[0]) });
        }
        for (const auto& i : indexBufferFirst)
        {
            indexBuffer.push_back(i + offsetFirst);
        }
 
        unsigned int offsetLast = static_cast<unsigned int>(vertexBuffer.size());
        for (const auto& v : tempLast)
        {
            vertexBuffer.push_back({ v.position, v.normalAndOpacity, wallValueToColor(controlPoints[controlPoints.size() - 1]) });
        }
        for (const auto& i : indexBufferLast)
        {
            indexBuffer.push_back(i + offsetLast);
        }
    }
}
 
static glm::mat4x4 shadowMatrix(glm::vec4 groundplane, glm::vec4 lightpos)
{
    // Find dot product between light position vector and ground plane normal.
    float dot = glm::dot(groundplane, lightpos);
 
    glm::mat4x4 shadowMat;
    shadowMat[0][0] = dot - lightpos.x * groundplane.x;
    shadowMat[1][0] = 0.f - lightpos.x * groundplane.y;
    shadowMat[2][0] = 0.f - lightpos.x * groundplane.z;
    shadowMat[3][0] = 0.f - lightpos.x * groundplane.w;
 
    shadowMat[0][1] = 0.f - lightpos.y * groundplane.x;
    shadowMat[1][1] = dot - lightpos.y * groundplane.y;
    shadowMat[2][1] = 0.f - lightpos.y * groundplane.z;
    shadowMat[3][1] = 0.f - lightpos.y * groundplane.w;
 
    shadowMat[0][2] = 0.f - lightpos.z * groundplane.x;
    shadowMat[1][2] = 0.f - lightpos.z * groundplane.y;
    shadowMat[2][2] = dot - lightpos.z * groundplane.z;
    shadowMat[3][2] = 0.f - lightpos.z * groundplane.w;
 
    shadowMat[0][3] = 0.f - lightpos.w * groundplane.x;
    shadowMat[1][3] = 0.f - lightpos.w * groundplane.y;
    shadowMat[2][3] = 0.f - lightpos.w * groundplane.z;
    shadowMat[3][3] = dot - lightpos.w * groundplane.w;
 
    return shadowMat;
}
 
void MultiphaseFlowComponent::generateMultiphaseFlowMesh(Mesh* mesh, const std::vector<ControlPoint>& controlPoints, MeshType liquidType)
{
    // Generate vertex and buffers for the main mesh.
    if (liquidType != MeshType::Wall && liquidType != MeshType::Shadow)
    {
        std::vector<Vertex> vertexBuffer;
        std::vector<unsigned int> indexBuffer;
        generateExtrusionVertices(vertexBuffer, indexBuffer, controlPoints, liquidType);
 
        // Set vertex and index buffers in the mesh object.
        if (!vertexBuffer.empty())
        {
            if (liquidType == MeshType::StratifiedOil)
            {
                calculateNormals(vertexBuffer, indexBuffer);
            }
 
            mesh->setVertexData(vertexBuffer.data(), vertexBuffer.size());
            mesh->setIndexes(std::move(indexBuffer));
        }
        else
        {
            mesh->clear();
        }
 
        // Update bounding box.
        for (auto& p : vertexBuffer)
        {
            geometryMin = glm::min(geometryMin, p.position);
            geometryMax = glm::max(geometryMax, p.position);
        }
    }
    else
    {
        std::vector<WallVertex> vertexBufferWall;
        std::vector<unsigned int> indexBuffer;
        generateExtrusionVerticesWall(vertexBufferWall, indexBuffer, controlPoints, liquidType);
 
        if (liquidType == MeshType::Shadow)
        {
            // Project shadow mesh on a plane.
            float l = glm::max(glm::max(glm::max(geometryMax.x, geometryMax.y), geometryMax.z), 1.f);
            glm::mat4x4 shadowMat = shadowMatrix(glm::vec4(0, 0, -1, geometryMin.z), glm::vec4(l* 1.5f, -l * 1.5f, l * 4.f, 1.f));
 
            for (auto& v : vertexBufferWall)
            {
                auto transformed = shadowMat * glm::vec4(v.position.x, v.position.y, v.position.z, 1);
                v.position = glm::vec3(transformed.x / transformed.w, transformed.y / transformed.w, transformed.z / transformed.w);
            }
        }
 
        // Set vertex and index buffers in the mesh object.
        if (!vertexBufferWall.empty())
        {
            mesh->setVertexData(vertexBufferWall.data(), vertexBufferWall.size());
            mesh->setIndexes(std::move(indexBuffer));
        }
        else
        {
            mesh->clear();
        }
 
 
        // Update bounding box.
        if (liquidType == MeshType::Wall)
        {
            for (auto& p : vertexBufferWall)
            {
                geometryMin = glm::min(geometryMin, p.position);
                geometryMax = glm::max(geometryMax, p.position);
            }
        }
    }
}
 
// Stretches the trajectory points to avoid overlapping of the mesh in the bends.
static void stretch(std::vector<MultiphaseFlowComponent::TrajectoryPoint>& points)
{
    assert(points.size() >= 2);
 
    glm::vec3 direction = points[1].pos - points[0].pos;
    glm::vec3 dirPrev = glm::normalize(direction);
    glm::vec3 prevPoint;
    glm::vec3 trajectoryOffset(0, 0, 0);
 
    for (size_t i = 0; i < points.size(); ++i)
    {
        // Calculate direction vector.
        if (i > 0)
        {
            if (i <= points.size() - 2)
            {
                direction = points[i + 1].pos - prevPoint;
            }
            else
            {
                direction = points[i].pos - prevPoint;
            }
        }
 
        float dirLength = glm::length(direction);
        if (dirLength > 0)
        {
            // Normalize direction.
            direction.x /= dirLength;
            direction.y /= dirLength;
            direction.z /= dirLength;
 
            // Save previous point.
            prevPoint = points[i].pos;
 
            // Get radius.
            const float r = points[i].radius;
 
            // Fix trajectory with sharp angles.
            // This will modify trajectory points and stretch the geometry at bends.
            auto dot = glm::dot(direction, dirPrev);
            dot = glm::clamp(glm::abs(dot), 0.f, 1.f);
            float angle = acos(dot);
            float ratio = angle / glm::half_pi<float>() * r * 0.65f;
 
            trajectoryOffset += (dirPrev + direction) * ratio;
            points[i].pos += trajectoryOffset;
 
            dirPrev = direction;
        }
    }
}
 
template <class T> class CatmullRomSpline
{
public:
    CatmullRomSpline(T x0, T x1, T x2, T x3, float t01, float t12, float t23, bool x0HasValue, bool x3HasValue)
    {
        T x0Scaled;
        if (x0HasValue)
        {
            auto delta = x1 - x0;
            delta /= t01;
            x0Scaled = x1 - delta * t12;
        }
        else
        {
            auto delta = x2 - x1;
            x0Scaled = x1 - delta;
        }
 
        T x3Scaled;
        if (x3HasValue)
        {
            auto delta = x2 - x3;
            delta /= t23;
            x3Scaled = x2 - delta * t12;
        }
        else
        {
            auto delta = x1 - x2;
            x3Scaled = x2 - delta;
        }
 
        // Compute control data
        _c0 = x1;
        _c1 = -x0Scaled * 0.5f + x2 * 0.5f;
        _c2 = x0Scaled - x1 * 2.5f + x2 * 2.0f - x3Scaled * 0.5f;
        _c3 = -x0Scaled * 0.5f + x1 * 1.5f - x2 * 1.5f + x3Scaled * 0.5f;
    }
 
    T Interpolate(float t)
    {
        return (((_c3 * t + _c2) * t + _c1) * t + _c0);
    }
 
private:
    T _c0;
    T _c1;
    T _c2;
    T _c3;
};
 
static void interpolate(std::vector<MultiphaseFlowComponent::TrajectoryPoint>& points)
{
    assert(points.size() >= 2);
 
    const float MAX_STEPS = glm::min(static_cast<float>(50000 / points.size()), 20.f);
    const float MIN_STEPS = glm::max(glm::floor(MAX_STEPS / 4), 1.f);
 
    if (MAX_STEPS <= 1)
    {
        return;
    }
 
    std::vector<MultiphaseFlowComponent::TrajectoryPoint> newPoints;
    size_t reservedSize = points.size() * static_cast<size_t>((MAX_STEPS + MIN_STEPS) / 2.f);
    newPoints.reserve(reservedSize);
 
    glm::vec3 direction = points[1].pos - points[0].pos;
    glm::vec3 dirPrev = glm::normalize(direction);
 
    for (size_t i = 0; i < points.size() - 1; i++)
    {
        if (i > 0)
        {
            if (i <= points.size() - 2)
            {
                direction = points[i + 1].pos - points[i - 1].pos;
            }
            else
            {
                direction = points[i].pos - points[i - 1].pos;
            }
        }
        direction = glm::normalize(direction);
 
        bool hasFirstPoint = (i > 0);
        bool hasLastPoint = (i < points.size() - 2);
 
        glm::vec3 p0 = hasFirstPoint ? points[i - 1].pos : glm::vec3();
        glm::vec3 f0 = hasFirstPoint ? points[i - 1].flowFraction : glm::vec3();
 
        const auto& p1 = points[i].pos;
        const auto& f1 = points[i].flowFraction;
 
        glm::vec3 p2 = points[i + 1].pos;
        glm::vec3 f2 = points[i + 1].flowFraction;
 
        glm::vec3 p3 = hasLastPoint ? points[i + 2].pos : glm::vec3();
        glm::vec3 f3 = hasLastPoint ? points[i + 2].flowFraction : glm::vec3();
 
        float t01 = hasFirstPoint ? glm::distance(p0, p1) : 0;
        float t12 = glm::distance(p1, p2);
        float t23 = hasLastPoint ? glm::distance(p2, p3) : 0;
 
        // Calculate number of interpolation steps based on the trajectory.
        auto dot = glm::dot(direction, dirPrev);
        auto angle = glm::acos(glm::clamp(glm::abs(dot), 0.f, 1.f));
        auto ratio = glm::clamp(angle / glm::half_pi<float>(), 0.f, 1.f);
 
        int steps = static_cast<int>(glm::lerp(MIN_STEPS, MAX_STEPS, ratio));
        assert(steps > 0);
 
        // Create spline objects.
        CatmullRomSpline<glm::vec3> splineTrajectory(p0, p1, p2, p3, t01, t12, t23, hasFirstPoint, hasLastPoint);
        CatmullRomSpline<glm::vec3> splineFlowFraction(f0, f1, f2, f3, t01, t12, t23, hasFirstPoint, hasLastPoint);
 
        // Interpolate.
        for (int j = 0; j < steps; j++)
        {
            float ratio_ = static_cast<float>(j) / static_cast<float>(steps);
 
            MultiphaseFlowComponent::TrajectoryPoint newPoint;
            newPoint.pos = splineTrajectory.Interpolate(ratio_);
            newPoint.radius = glm::lerp(points[i].radius, points[i + 1].radius, ratio_);
            newPoint.flowFraction = splineFlowFraction.Interpolate(ratio_);
            newPoint.wallValue = glm::lerp(points[i].wallValue, points[i + 1].wallValue, ratio_);
            newPoint.md = glm::lerp(points[i].md, points[i + 1].md, ratio_);
 
            newPoints.push_back(newPoint);
        }
 
        dirPrev = direction;
    }
 
    // Add last point.
    newPoints.push_back(points[points.size() - 1]);
 
    // Swap data.
    std::swap(points, newPoints);
}
 
void MultiphaseFlowComponent::updateGroundGrid()
{
    // Set visibility.
    if (groundGridMeshEntity)
    {
        auto meshSceneCompoment = groundGridMeshEntity->getComponent<MeshRenderComponent>();
        meshSceneCompoment->setVisible(showGrid);
    }
 
    // Exit if grid is not visible.
    if (!showGrid)
    {
        return;
    }
 
    glm::vec3 center = (geometryMax + geometryMin) * 0.5f;
    float edgeLength = glm::max(geometryMax.x - geometryMin.x, geometryMax.y - geometryMin.y) * 1.5f;           // Make the grid 150% bigger than the actual mesh.
    glm::vec3 corner = center - edgeLength * 0.5f;
 
    // Create grid
    std::vector<PositionVertex> vertexBufferGrid;
    std::vector<unsigned int> indexBufferGrid;
    static const int NUM_OF_SECTORS = 10;
    for (int j = 0; j < 2; j++)
    {
        for (int i = 0; i <= NUM_OF_SECTORS; i++)
        {
            PositionVertex v0;
            PositionVertex v1;
 
            if (j == 0)
            {
                v0.position = glm::vec3(corner.x + edgeLength / static_cast<float>(NUM_OF_SECTORS)* static_cast<float>(i), corner.y, geometryMin.z * 1.01f);
                v1.position = glm::vec3(corner.x + edgeLength / static_cast<float>(NUM_OF_SECTORS)* static_cast<float>(i), corner.y + edgeLength, geometryMin.z * 1.01f);
            }
            else
            {
                v0.position = glm::vec3(corner.x, corner.y + edgeLength / static_cast<float>(NUM_OF_SECTORS)* static_cast<float>(i), geometryMin.z * 1.01f);
                v1.position = glm::vec3(corner.x + edgeLength, corner.y + edgeLength / static_cast<float>(NUM_OF_SECTORS)* static_cast<float>(i), geometryMin.z * 1.01f);
            }
 
            vertexBufferGrid.push_back(v0);
            vertexBufferGrid.push_back(v1);
 
            indexBufferGrid.push_back(static_cast<unsigned int>(vertexBufferGrid.size()) - 2);
            indexBufferGrid.push_back(static_cast<unsigned int>(vertexBufferGrid.size()) - 1);
        }
    }
 
    // Update mesh with the new vertex and index buffers.
    Mesh* meshGrid = context->meshManager->get(groundGridMeshEntity->getComponent<MeshComponent>()->meshHandle);
    meshGrid->setVertexData(vertexBufferGrid.data(), vertexBufferGrid.size());
    meshGrid->setIndexes(std::move(indexBufferGrid));
    meshGrid->primitiveType = Cogs::PrimitiveType::LineList;
}
 
void MultiphaseFlowComponent::updateIndicatorRing(const std::vector<TrajectoryPoint>& trajectory)
{
    // Set visibility.
    if (indicatorRingEntity)
    {
        auto meshSceneCompoment = indicatorRingEntity->getComponent<MeshRenderComponent>();
        meshSceneCompoment->setVisible(showIndicatorRing);
    }
 
    // Exit if the ring is not visible.
    if (!showIndicatorRing)
    {
        return;
    }
 
    // Find position.
    assert(trajectory.size() > 1);
    glm::vec3 positionPoint;
    float radius = 0.f;
    glm::quat rotation;
    for (size_t i = 1; i < trajectory.size(); i++)
    {
        if (indicatorRingPosition < trajectory[i].md || i+1 == trajectory.size())
        {
            auto ratio = (trajectory[i].md - indicatorRingPosition) / (trajectory[i].md - trajectory[i - 1].md);
            ratio = glm::clamp(ratio, 0.f, 1.f);
 
            positionPoint = glm::lerp(trajectory[i].pos, trajectory[i - 1].pos, ratio);
            radius = glm::lerp(trajectory[i].radius, trajectory[i - 1].radius, ratio) * 1.3f;
            rotation = glm::slerp(trajectory[i].rotation, trajectory[i - 1].rotation, ratio);
 
            break;
        }
    }
 
    std::vector<PositionVertex> vertexBuffer;
    std::vector<unsigned int> indexBuffer;
    static const int NUM_OF_STEPS = 32;
    auto angle = 0.0;
    auto angleStep = glm::two_pi<double>() / static_cast<double>(NUM_OF_STEPS);
 
    for (auto i = 0; i < NUM_OF_STEPS; i++)
    {
        auto s = static_cast<float>(sin(angle)) * radius;
        auto c = static_cast<float>(cos(angle)) * radius;
        angle += angleStep;
 
        PositionVertex v;
        v.position = rotation * glm::vec3(s, c, 0);
        v.position += positionPoint;
 
        vertexBuffer.push_back(v);
 
        if (i == NUM_OF_STEPS - 1)
        {
            indexBuffer.push_back(i);
            indexBuffer.push_back(0);
        }
        else
        {
            indexBuffer.push_back(i);
            indexBuffer.push_back(i + 1);
        }
    }
 
    // Update mesh with the new vertex and index buffers.
    Mesh* meshGrid = context->meshManager->get(indicatorRingEntity->getComponent<MeshComponent>()->meshHandle);
    meshGrid->setVertexData(vertexBuffer.data(), vertexBuffer.size());
    meshGrid->setIndexes(std::move(indexBuffer));
    meshGrid->primitiveType = Cogs::PrimitiveType::LineList;
 
    // Update material color.
    auto meshRenderComp = indicatorRingEntity->getComponent<MeshRenderComponent>();
    meshRenderComp->material->setVec4Property(DefaultMaterial::DiffuseColor, glm::vec4(indicatorRingColor, 1));
}
 
void MultiphaseFlowComponent::updateAxisAnnotations()
{
    Mesh* mesh = context->meshManager->get(annotationAxisMeshEntity->getComponent<MeshComponent>()->meshHandle);
 
    // If axial extents are hidden.
    if (!showAxialExtents)
    {
        // Hide axial extents meshes.
        mesh->clear();
 
        // Hide text.
        for (int i = AnnotationAxis::AxisFirst; i <= AnnotationAxis::AxisLast; i++)
        {
            auto textComp = textEntity[i]->getComponent<TextComponent>();
            textComp->texts.clear();
        }
 
        // Exit.
        return;
    }
 
    glm::vec3 center = (geometryMax + geometryMin) * 0.5f;
    float edgeLength = glm::max(geometryMax.x - geometryMin.x, geometryMax.y - geometryMin.y) * 1.5f;           // Make the grid 150% bigger than the actual mesh.
 
    std::vector<PositionVertex> vertexBuffer;
    std::vector<unsigned int> indexBuffer;
 
    // Calculate trajectory bounding box.
    glm::vec3 trajectoryMin = trajectoryPoints[0];
    glm::vec3 trajectoryMax = trajectoryPoints[0];
    for (const auto& p : trajectoryPoints)
    {
        trajectoryMin = glm::min(trajectoryMin, p);
        trajectoryMax = glm::max(trajectoryMax, p);
    }
 
    // Show annotations based on the trajectory bounding box.
    bool xShowAnnotation = glm::abs(trajectoryMax.x - trajectoryMin.x) > 0.1f;
    bool yShowAnnotation = glm::abs(trajectoryMax.y - trajectoryMin.y) > 0.1f;
    bool zShowAnnotation = glm::abs(trajectoryMax.z - trajectoryMin.z) > 0.1f;
 
    // X line.
    float xLinePosZ = geometryMax.z + edgeLength * 0.01f;
    if (xShowAnnotation)
    {
        vertexBuffer.push_back({ glm::vec3(scaledTrajectoryMin.x, center.y, xLinePosZ) });
        vertexBuffer.push_back({ glm::vec3(scaledTrajectoryMax.x, center.y, xLinePosZ) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
        vertexBuffer.push_back({ glm::vec3(scaledTrajectoryMin.x, center.y, xLinePosZ - edgeLength * 0.005f) });
        vertexBuffer.push_back({ glm::vec3(scaledTrajectoryMin.x, center.y, xLinePosZ + edgeLength * 0.005f) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
        vertexBuffer.push_back({ glm::vec3(scaledTrajectoryMax.x, center.y, xLinePosZ - edgeLength * 0.005f) });
        vertexBuffer.push_back({ glm::vec3(scaledTrajectoryMax.x, center.y, xLinePosZ + edgeLength * 0.005f) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
    }
 
    // Y line.
    float yLinePosX = geometryMin.x - edgeLength * 0.01f;
    if (yShowAnnotation)
    {
        vertexBuffer.push_back({ glm::vec3(yLinePosX, scaledTrajectoryMin.y, center.z) });
        vertexBuffer.push_back({ glm::vec3(yLinePosX, scaledTrajectoryMax.y, center.z) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
        vertexBuffer.push_back({ glm::vec3(yLinePosX - edgeLength * 0.005f, scaledTrajectoryMin.y, center.z) });
        vertexBuffer.push_back({ glm::vec3(yLinePosX + edgeLength * 0.005f, scaledTrajectoryMin.y, center.z) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
        vertexBuffer.push_back({ glm::vec3(yLinePosX - edgeLength * 0.005f, scaledTrajectoryMax.y, center.z) });
        vertexBuffer.push_back({ glm::vec3(yLinePosX + edgeLength * 0.005f, scaledTrajectoryMax.y, center.z) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
    }
 
    // Z line.
    float zLinePosX = geometryMax.x + edgeLength * 0.01f;
    if (zShowAnnotation)
    {
        vertexBuffer.push_back({ glm::vec3(zLinePosX, center.y, scaledTrajectoryMin.z) });
        vertexBuffer.push_back({ glm::vec3(zLinePosX, center.y, scaledTrajectoryMax.z) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
        vertexBuffer.push_back({ glm::vec3(zLinePosX - edgeLength * 0.005f, center.y, scaledTrajectoryMin.z) });
        vertexBuffer.push_back({ glm::vec3(zLinePosX + edgeLength * 0.005f, center.y, scaledTrajectoryMin.z) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
        vertexBuffer.push_back({ glm::vec3(zLinePosX - edgeLength * 0.005f, center.y, scaledTrajectoryMax.z) });
        vertexBuffer.push_back({ glm::vec3(zLinePosX + edgeLength * 0.005f, center.y, scaledTrajectoryMax.z) });
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 2);
        indexBuffer.push_back(static_cast<unsigned int>(vertexBuffer.size()) - 1);
    }
 
    // Update mesh with the new vertex and index buffers.
    mesh->setVertexData(vertexBuffer.data(), vertexBuffer.size());
    mesh->setIndexes(std::move(indexBuffer));
    mesh->primitiveType = Cogs::PrimitiveType::LineList;
 
    // Clear text.
    TextComponent* textComp[AnnotationAxis::AxisLast + 1] = { nullptr, nullptr, nullptr };
    for (int i = AnnotationAxis::AxisFirst; i <= AnnotationAxis::AxisLast; i++)
    {
        textComp[i] = textEntity[i]->getComponent<TextComponent>();
        textComp[i]->color = glm::vec4(0.f, 0.f, 0.f, 1.f);
        textComp[i]->texts.clear();
        textComp[i]->positions.clear();
        textComp[i]->positionMode = PositionMode::World;
        textComp[i]->setChanged();
        textComp[i]->getComponent<SceneComponent>()->visible = true;
    }
 
    // X axis annotation.
    if (xShowAnnotation)
    {
        std::stringstream xMinText;
        xMinText << static_cast<int>(floor(trajectoryMin.x * unitScale)) << " " << unitText;
        textComp[AnnotationAxis::X]->texts.push_back(xMinText.str());
        textComp[AnnotationAxis::X]->positions.push_back(glm::vec3(scaledTrajectoryMin.x, center.y, xLinePosZ + edgeLength * 0.01f));
 
        std::stringstream xMaxText;
        xMaxText << static_cast<int>(floor(trajectoryMax.x * unitScale)) << " " << unitText;
        textComp[AnnotationAxis::X]->texts.push_back(xMaxText.str());
        textComp[AnnotationAxis::X]->positions.push_back(glm::vec3(scaledTrajectoryMax.x, center.y, xLinePosZ + edgeLength * 0.01f));
 
        textComp[AnnotationAxis::X]->texts.push_back("X");
        textComp[AnnotationAxis::X]->positions.push_back(glm::vec3((scaledTrajectoryMax.x + scaledTrajectoryMin.x) / 2, center.y, xLinePosZ + edgeLength * 0.005f));
    }
 
    // Y axis annotation.
    if (yShowAnnotation)
    {
        std::stringstream yMinText;
        yMinText << static_cast<int>(floor(trajectoryMin.y * unitScale)) << " " << unitText;
        textComp[AnnotationAxis::YMin]->texts.push_back(yMinText.str());
        textComp[AnnotationAxis::YMin]->positions.push_back(glm::vec3(yLinePosX - edgeLength * 0.01f, scaledTrajectoryMin.y, center.z));
 
        std::stringstream yMaxText;
        yMaxText << static_cast<int>(floor(trajectoryMax.y * unitScale)) << " " << unitText;
        textComp[AnnotationAxis::YMax]->texts.push_back(yMaxText.str());
        textComp[AnnotationAxis::YMax]->positions.push_back(glm::vec3(yLinePosX - edgeLength * 0.01f, scaledTrajectoryMax.y, center.z));
 
        textComp[AnnotationAxis::Y]->texts.push_back("Y");
        textComp[AnnotationAxis::Y]->positions.push_back(glm::vec3(yLinePosX + edgeLength * 0.005f, (scaledTrajectoryMin.y + scaledTrajectoryMax.y) / 2, center.z));
    }
 
    // Z axis annotation.
    if (zShowAnnotation)
    {
        std::stringstream zMaxText;
        zMaxText << "  " << static_cast<int>(ceil(trajectoryMax.z * unitScale)) << " " << unitText << "  ";
        textComp[AnnotationAxis::Z]->texts.push_back(zMaxText.str());
        textComp[AnnotationAxis::Z]->positions.push_back(glm::vec3(zLinePosX + edgeLength * 0.005f, center.y, scaledTrajectoryMax.z));
 
        std::stringstream zMinText;
        zMinText << "  " << static_cast<int>(floor(trajectoryMin.z * unitScale)) << " " << unitText << "  ";
        textComp[AnnotationAxis::Z]->texts.push_back(zMinText.str());
        textComp[AnnotationAxis::Z]->positions.push_back(glm::vec3(zLinePosX + edgeLength * 0.005f, center.y, scaledTrajectoryMin.z));
 
        textComp[AnnotationAxis::Z]->texts.push_back("Z");
        textComp[AnnotationAxis::Z]->positions.push_back(glm::vec3(zLinePosX + edgeLength * 0.005f, center.y, (scaledTrajectoryMin.z + scaledTrajectoryMax.z) / 2));
    }
 
    // Text has been changed.
    for (int i = AnnotationAxis::AxisFirst; i <= AnnotationAxis::AxisLast; i++)
    {
        textComp[i]->setChanged();
    }
}
 
void MultiphaseFlowComponent::removeDuplicatedPoints(std::vector<TrajectoryPoint>& points)
{
    assert(!points.empty());
    std::vector<TrajectoryPoint> newPoints;
    newPoints.reserve(points.size());
    newPoints.push_back(points[0]);
 
    for (size_t i = 1; i < points.size(); i++)
    {
        if (points[i].pos != newPoints[newPoints.size() - 1].pos)
        {
            newPoints.push_back(points[i]);
        }
    }
 
    points.swap(newPoints);
}
 
void MultiphaseFlowComponent::applyScaling(std::vector<TrajectoryPoint>& points)
{
    assert(!points.empty());
 
    float bottomHeight = points[0].pos.z;
    for (auto& p : points)
    {
        bottomHeight = glm::min(bottomHeight, p.pos.z);
    }
 
    for (auto& p : points)
    {
        p.pos.z = (p.pos.z - bottomHeight) * heightScale;
    }
 
    for (auto& p : points)
    {
        p.radius *= radiusScale;
    }
}
 
void MultiphaseFlowComponent::center(std::vector<TrajectoryPoint>& points)
{
    assert(!points.empty());
 
    glm::vec3 min = points[0].pos;
    glm::vec3 max = points[0].pos;
    for (auto& p : points)
    {
        min = glm::min(min, p.pos);
        max = glm::max(max, p.pos);
    }
 
    glm::vec3 center;
    center.x = (min.x + max.x) / 2;
    center.y = (min.y + max.y) / 2;
    center.z = (min.z + max.z) / 2;
 
    for (auto& p : points)
    {
        p.pos -= center;
    }
}
 
static float crosssectionAreaToHeight(float alpha)
{
    static const float fac = glm::pow(1.5f * glm::pi<float>(), 1.f / 3.f);
    const float beta = glm::pi<float>() * alpha + fac * (1 - 2 * alpha + glm::pow(alpha, 1.f / 3.f) - glm::pow(1 - alpha, 1.f / 3.f));
    float res = 0.5f - 0.5f * glm::cos(beta);
    return res;
}
 
static glm::quat myRotation(glm::vec3 dir, glm::vec3 up = glm::vec3(0, 1, 0))
{
    glm::vec3 x = glm::cross(up, dir);
    x = glm::normalize(x);
 
    glm::vec3 y = glm::cross(dir, x);
    y = glm::normalize(y);
 
    glm::mat3x3 cMatrix(x, y, dir);
 
    return glm::quat(cMatrix);
}
 
void MultiphaseFlowComponent::finalizeControlPoints(std::vector<ControlPoint>& controlPoints)
{
    static const int DISTRIBUTION_SIZE = 5;
    static float distributionTable[DISTRIBUTION_SIZE + 1];
    static bool distributionTableInitialized = false;
    if (!distributionTableInitialized)
    {
        for (int x = 0; x <= DISTRIBUTION_SIZE; x++)
        {
            float y = expf(-(x * x * 0.1f));
            distributionTable[x] = y;
        }
        distributionTableInitialized = true;
    }
 
    for (int i = 0; i < static_cast<int>(controlPoints.size()); i++)
    {
        auto& controlPoint = controlPoints[i];
 
        // Interpolate annular flow.
        controlPoint.annularFlowRatio = 0.f;
 
        float sum = 0.f;
        for (int j = -DISTRIBUTION_SIZE; j <= DISTRIBUTION_SIZE; j++)
        {
            int k = i + j;
            if (k < 0) k = 0;
            if (k > static_cast<int>(controlPoints.size() - 1)) k = static_cast<int>(controlPoints.size() - 1);
 
            auto& sibling = controlPoints[k];
 
            float ratio = distributionTable[j < 0 ? -j : j];
            controlPoint.annularFlowRatio += sibling.annularFlow ? ratio : 0.f;
            sum += ratio;
        }
 
        assert(sum > 0);
        controlPoint.annularFlowRatio /= sum;
        controlPoint.annularFlowRatio = glm::clamp(controlPoint.annularFlowRatio, 0.f, 1.f);
 
        // Set binary inverted flow.
        controlPoint.invertedFlowRatio = controlPoint.invertedFlow ? 1.f : 0.f;
    }
 
    // Get ranges for annular flow.
    std::vector<std::pair<int, int>> annularRanges;
    int start = -1;
    for (int i = 0; i < static_cast<int>(controlPoints.size()); i++)
    {
        auto& controlPoint = controlPoints[i];
        if (start == -1 && controlPoint.annularFlowRatio > 0.f)
        {
            start = glm::max(i - 1, 0);
        }
        else
            if (start >= 0 && controlPoint.annularFlowRatio == 0.f)
            {
                assert(i > 0);
                annularRanges.push_back(std::pair<int, int>(start, i));
                start = -1;
            }
    }
    if (start >= 0)
    {
        annularRanges.push_back(std::pair<int, int>(start, static_cast<int>(controlPoints.size() - 1)));
    }
 
    // Interpolate rotations and inverted flow.
    for (int i = 0; i < (int)annularRanges.size(); i++)
    {
        auto& p0 = controlPoints[annularRanges[i].first];
        auto& p1 = controlPoints[annularRanges[i].second];
 
        float p0InvertedFlow = p0.invertedFlow ? 1.f : 0.f;
        float p1InvertedFlow = p1.invertedFlow ? 1.f : 0.f;
 
        int dist = annularRanges[i].second - annularRanges[i].first;
        if (dist > 1)
        {
            for (int j = 0; j <= dist; j++)
            {
                float ratio = static_cast<float>(j) / static_cast<float>(dist);
                ratio = sCurve(ratio);
                ratio = glm::clamp(ratio, 0.f, 1.f);
 
                auto up = glm::lerp(p0.rotationUpVector, p1.rotationUpVector, ratio);
 
                int k = annularRanges[i].first + j;
                controlPoints[k].rotation = myRotation(controlPoints[k].direction, up);
 
                // Interpolate inverted flow.
                controlPoints[k].invertedFlowRatio = glm::lerp(p0InvertedFlow, p1InvertedFlow, ratio);
            }
        }
    }
}
 
void MultiphaseFlowComponent::updateShadowVisibility()
{
    auto* cam = context->cameraSystem->getMainCamera();
    auto transform = static_cast<TransformComponent*>(cam->getComponent<TransformComponent>());
 
    bool show = (transform->position.z > geometryMin.z) && (trajectoryPoints.size() == radius.size()) && showShadow;
 
    if (meshEntity[MeshType::Shadow])
    {
        auto meshSceneCompoment = meshEntity[MeshType::Shadow]->getComponent<MeshRenderComponent>();
        if (meshSceneCompoment)
        {
            meshSceneCompoment->setVisible(show);
        }
    }
}
 
void MultiphaseFlowComponent::updateMeshes()
{
    assert(updateNeeded);
 
    updateNeeded = false;
 
    assert(flowFraction.empty() || flowFraction.size() == trajectoryPoints.size());
    assert(wallValues.empty() || wallValues.size() == trajectoryPoints.size());
 
    std::vector<TrajectoryPoint> trajectory(trajectoryPoints.size());
    float md = 0.0f;
    for (size_t i = 0; i < trajectoryPoints.size(); i++)
    {
        auto& trajPoint = trajectory[i];
        trajPoint.pos = trajectoryPoints[i];
        trajPoint.radius = radius[i];
        trajPoint.flowFraction = !flowFraction.empty() ? flowFraction[i] : glm::vec3();
        trajPoint.wallValue = !wallValues.empty() ? wallValues[i] : 0.0f;
 
        if (i > 0)
        {
            md += glm::distance(trajectory[i - 1].pos, trajectory[i].pos);
        }
 
        trajPoint.md = md;
    }
 
    removeDuplicatedPoints(trajectory);
 
    if (centerMesh)
    {
        center(trajectory);
    }
 
    applyScaling(trajectory);
 
    stretch(trajectory);
    interpolate(trajectory);
    stretch(trajectory);
 
    if (centerMesh)
    {
        center(trajectory);
    }
    else
    {
        auto posOffset = trajectoryPoints[0] - trajectory[0].pos;
        for (auto &v : trajectory)
        {
            v.pos += posOffset;
        }
    }
 
    std::vector<ControlPoint> controlPoints;
    controlPoints.resize(trajectory.size());
    bool visualizeMultiphaseFlow = !flowFraction.empty();
    scaledTrajectoryMin = trajectory[0].pos;
    scaledTrajectoryMax = trajectory[0].pos;
    bool waterInMultiphaseFlow = false;
    bool oilInMultiphaseFlow = false;
    bool gasInMultiphaseFlow = false;
    for (size_t i = 0; i < trajectory.size(); i++)
    {
        auto& controlPoint = controlPoints[i];
        controlPoint.position = trajectory[i].pos;
        controlPoint.radius = trajectory[i].radius;
        controlPoint.wall = glm::clamp(!wallValues.empty() ? trajectory[i].wallValue : 0.f, 0.f, 1.f);
 
        // Update trajectory AABB
        scaledTrajectoryMin = glm::min(scaledTrajectoryMin, trajectory[i].pos);
        scaledTrajectoryMax = glm::max(scaledTrajectoryMax, trajectory[i].pos);
 
        // Calculate direction vector.
        controlPoint.direction = trajectory[i == 0 ? 0 : i - 1].pos - trajectory[i < trajectory.size() - 1 ? i + 1 : i].pos;
        controlPoint.direction = glm::normalize(controlPoint.direction);
 
        // Segment rotation
        auto d = glm::dot(controlPoint.direction, glm::vec3(0, 1, 0));
        auto right = glm::abs(d) > 0.95f ? glm::vec3(1, 0, 0) : glm::vec3(0, 1, 0);
        if (i == 0)
        {
            controlPoint.rotationUpVector = right;
        }
        else
        {
            auto rotPos = myRotation(controlPoint.direction, right);
            auto rotNeg = myRotation(controlPoint.direction, -right);
 
            auto upPos = rotPos * glm::vec3(1, 0, 0);
            auto upNeg = rotNeg * glm::vec3(1, 0, 0);
 
            auto upPrev = controlPoints[i - 1].rotation * glm::vec3(1, 0, 0);
 
            auto dotPos = glm::dot(upPos, upPrev);
            auto dotNeg = glm::dot(upNeg, upPrev);
 
            controlPoint.rotationUpVector = right * (dotPos > dotNeg ? 1.0f : -1.0f);
        }
        controlPoint.rotation = myRotation(controlPoint.direction, controlPoint.rotationUpVector);
        trajectory[i].rotation = controlPoint.rotation;
 
        // Pre-compute data used to show liquid fraction.
        auto& flowData = controlPoint.flow;
        flowData.water = visualizeMultiphaseFlow ? glm::max(trajectory[i].flowFraction.x, 0.f) : 0.f;
        flowData.oil = visualizeMultiphaseFlow ? glm::max(trajectory[i].flowFraction.y, 0.f) : 0.f;
        flowData.gas = visualizeMultiphaseFlow ? glm::max(trajectory[i].flowFraction.z, 0.f) : 0.f;
 
        // Re-normalize amount of liquids and gas.
        float total = flowData.water + flowData.oil + flowData.gas;
        if (total > 0.f)
        {
            flowData.water /= total;
            flowData.oil /= total;
            flowData.gas /= total;
        }
 
        // Check visibility.
        if (flowData.water > 0.f) { waterInMultiphaseFlow = true; }
        if (flowData.oil > 0.f) { oilInMultiphaseFlow = true; }
        if (flowData.gas > 0.f) { gasInMultiphaseFlow = true; }
 
        // Calculate liquid height in the pipe based on the percentage of the pipe area.
        float h1 = crosssectionAreaToHeight(flowData.water);
        float h2 = crosssectionAreaToHeight(flowData.water + flowData.oil);
 
        flowData.heightWaterOil = -1 * (h1 - 0.5f) * 2.f;
        flowData.heightOilGas = -1 * (h2 - 0.5f) * 2.f;
 
        flowData.widthWaterOil = glm::sqrt(1 - (1 - h1 * 2) * (1 - h1 * 2));
        flowData.widthOilGas = glm::sqrt(1 - (1 - h2 * 2) * (1 - h2 * 2));
 
        // Check for inverted crosssection (water always at the bottom).
        glm::vec3 testVec = controlPoint.rotation * glm::vec3(1, 0, 0);
        controlPoint.invertedFlow = (testVec.z > 0.f);
 
        // Check for annular flow.
        bool isAnnularFlow = (std::fabs(controlPoint.direction.z) >= ANNULAR_FLOW_INCLINATION);
        controlPoint.annularFlow = isAnnularFlow;
    }
 
    // Interpolate control data.
    finalizeControlPoints(controlPoints);
 
    // Reset geometry AABB.
    geometryMin = glm::vec3(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
    geometryMax = glm::vec3(std::numeric_limits<float>::lowest(), std::numeric_limits<float>::lowest(), std::numeric_limits<float>::lowest());
 
    // Update all meshes.
    for (int i = MeshType::First; i <= MeshType::Last; i++)
    {
        auto liquidType = MeshType(i);
 
        // Get mesh.
        Mesh* mesh = context->meshManager->get(meshEntity[i]->getComponent<MeshComponent>()->meshHandle);
 
        // Update material.
        updateMaterial(meshEntity[i], liquidType);
 
        bool visible = true;
        switch (liquidType)
        {
        case MeshType::StratifiedWater:
        case MeshType::AnnularWater:
            visible = showWater && waterInMultiphaseFlow && visualizeMultiphaseFlow;
            break;
        case MeshType::StratifiedOil:
        case MeshType::AnnularOil:
            visible = showOil && oilInMultiphaseFlow && visualizeMultiphaseFlow;
            break;
        case MeshType::StratifiedGas:
        case MeshType::AnnularGas:
            visible = showGas && gasInMultiphaseFlow && visualizeMultiphaseFlow;
            break;
        case MeshType::Wall:
            visible = (wallValues.size() == trajectoryPoints.size() || !visualizeMultiphaseFlow);
            if (visible)
            {
                // Modify material based on the visibility of multiphase flow mesh.
                auto meshRenderComp = meshEntity[i]->getComponent<MeshRenderComponent>();
                auto material = meshRenderComp->material;
                material->setOption("CullMode", visualizeMultiphaseFlow ? "Back" : "Front");
 
                material->setVariant("EnableLighting", !visualizeMultiphaseFlow);
                material->setVariant("LightModel", visualizeMultiphaseFlow ? "BaseColor" : "Phong");
            }
            break;
        case MeshType::Shadow:
        {
            visible = (trajectoryPoints.size() == radius.size());
        }
        break;
        default:
            assert(false && "Unknown liquid type");
            break;
        }
 
        if (visible)
        {
            generateMultiphaseFlowMesh(mesh, controlPoints, liquidType);
        }
        else
        {
            mesh->clear();
        }
    }
 
    // Update grid, indicator ring and text.
    updateGroundGrid();
    updateAxisAnnotations();
    updateIndicatorRing(trajectory);
}
 
void MultiphaseFlowComponent::update()
{
    if (trajectoryPoints.size() < 2 || trajectoryPoints.size() != radius.size())
    {
        return;
    }
 
    if (showAxialExtents != currentShowAxialExtents ||
        showShadow != currentShowGrid ||
        showGrid != currentShowGrid ||
        centerMesh != currentCenterMesh ||
        showIndicatorRing != currentShowIndicatorRing)
    {
        currentShowAxialExtents = showAxialExtents;
        currentShowGrid = showShadow;
        currentShowGrid = showGrid;
        currentCenterMesh = centerMesh;
        currentShowIndicatorRing = showIndicatorRing;
        updateNeeded = true;
    }
 
    if (updateNeeded)
    {
        updateFontSize();
        updateMeshes();
    }
 
    updateShadowVisibility();
 
    // Set text horizontal alignment for Y an Z axis based on camera position.
    // Avoid text being displayed under the 3D model of a pipe.
    auto* cam = context->cameraSystem->getMainCamera();
    auto transform = static_cast<TransformComponent*>(cam->getComponent<TransformComponent>());
    auto* zTextComp = textEntity[AnnotationAxis::Z]->getComponent<TextComponent>();
    zTextComp->horizontalJustification = (transform->position.y <= 0 ? HorizontalJustification::Left : HorizontalJustification::Right);
 
    auto* yMinTextComp = textEntity[AnnotationAxis::YMin]->getComponent<TextComponent>();
    auto* yMaxTextComp = textEntity[AnnotationAxis::YMax]->getComponent<TextComponent>();
    auto* yTextComp = textEntity[AnnotationAxis::Y]->getComponent<TextComponent>();
 
    if (yMinTextComp != nullptr && yMaxTextComp != nullptr && yTextComp != nullptr)
    {
        yMinTextComp->horizontalJustification = (transform->position.y > 0 ? HorizontalJustification::Left : HorizontalJustification::Right);
        yMaxTextComp->horizontalJustification = (transform->position.y > 0 ? HorizontalJustification::Left : HorizontalJustification::Right);
        yTextComp->horizontalJustification = (transform->position.y <= 0 ? HorizontalJustification::Left : HorizontalJustification::Right);
    }
}
 
void MultiphaseFlowComponent::calculateNormals(std::vector<Vertex>& vertexBuffer, std::vector<unsigned int>& indexBuffer)
{
    for (size_t index = 0; index < indexBuffer.size() / 3; index++)
    {
        auto i0 = indexBuffer[index * 3 + 0];
        auto i1 = indexBuffer[index * 3 + 1];
        auto i2 = indexBuffer[index * 3 + 2];
 
        auto& p0 = vertexBuffer[i0];
        auto& p1 = vertexBuffer[i1];
        auto& p2 = vertexBuffer[i2];
 
        auto p = p1.position - p0.position;
        auto q = p2.position - p0.position;
        auto s = p2.position - p1.position;
 
        if (glm::length(p) > 0.0001f &&
            glm::length(q) > 0.0001f &&
            glm::length(s) > 0.0001f)
        {
            auto n = glm::cross(p, q);
            n = glm::normalize(n);
 
            p0.normalAndOpacity.x = n.x;
            p0.normalAndOpacity.y = n.y;
            p0.normalAndOpacity.z = n.z;
 
            p1.normalAndOpacity.x = n.x;
            p1.normalAndOpacity.y = n.y;
            p1.normalAndOpacity.z = n.z;
 
            p2.normalAndOpacity.x = n.x;
            p2.normalAndOpacity.y = n.y;
            p2.normalAndOpacity.z = n.z;
        }
    }
}