From 24d38d0296f8913b49e2b44975226906f0bdd7fa Mon Sep 17 00:00:00 2001 From: Hans-Christian Helltegen Date: Fri, 27 Jun 2014 15:18:34 -0400 Subject: [PATCH] Implemented tracing for Lorentz force trajectories --- .../rendering/renderablefieldlines.h | 6 +- include/openspace/util/kameleonwrapper.h | 17 +- src/rendering/renderablefieldlines.cpp | 87 ++++++---- src/util/kameleonwrapper.cpp | 158 +++++++++++++++++- 4 files changed, 220 insertions(+), 48 deletions(-) diff --git a/include/openspace/rendering/renderablefieldlines.h b/include/openspace/rendering/renderablefieldlines.h index e26c7c2b54..3f663c4f73 100644 --- a/include/openspace/rendering/renderablefieldlines.h +++ b/include/openspace/rendering/renderablefieldlines.h @@ -58,7 +58,7 @@ private: std::vector _filenames; std::vector _seedPoints; - ghoul::opengl::ProgramObject *_fieldlinesProgram, *_seedpointsProgram; + ghoul::opengl::ProgramObject *_fieldlinesProgram; GLuint _VAO, _seedpointVAO; std::mutex* _shaderMutex; @@ -66,8 +66,8 @@ private: ghoul::filesystem::File* _vertexSourceFile; ghoul::filesystem::File* _fragmentSourceFile; - std::vector _lineStart, _seedpointStart; - std::vector _lineCount, _seedpointCount; + std::vector _lineStart; + std::vector _lineCount; bool _programUpdateOnSave; void safeShaderCompilation(); diff --git a/include/openspace/util/kameleonwrapper.h b/include/openspace/util/kameleonwrapper.h index 4141bdec52..6b85680231 100644 --- a/include/openspace/util/kameleonwrapper.h +++ b/include/openspace/util/kameleonwrapper.h @@ -25,7 +25,6 @@ #ifndef KAMELEONWRAPPER_H_ #define KAMELEONWRAPPER_H_ -//#include #include namespace ccmc { @@ -43,12 +42,12 @@ public: BATSRUS // Magnetosphere }; - enum TraceDirection { + enum class TraceDirection { FORWARD = 1, BACK = -1 }; - enum FieldlineEnd { + enum class FieldlineEnd { NORTH, SOUTH, OUT @@ -60,15 +59,25 @@ public: float* getUniformSampledVectorValues(const std::string& xVar, const std::string& yVar, const std::string& zVar, glm::size3_t outDimensions); - std::vector > getFieldLines(const std::string& xVar, + std::vector > getClassifiedFieldLines(const std::string& xVar, const std::string& yVar, const std::string& zVar, std::vector seedPoints, float stepSize); + std::vector > getFieldLines(const std::string& xVar, + const std::string& yVar, const std::string& zVar, + std::vector seedPoints, float stepSize, glm::vec3 color); + + std::vector > getLorentzTrajectories(std::vector seedPoints, + glm::vec3 color, float stepsize); + private: std::vector traceCartesianFieldline(const std::string& xVar, const std::string& yVar, const std::string& zVar, glm::vec3 seedPoint, float stepSize, TraceDirection direction, FieldlineEnd& end); + std::vector traceLorentzTrajectory(glm::vec3 seedPoint, + float stepsize, float eCharge); + void getGridVariables(std::string& x, std::string& y, std::string& z); void progressBar(int current, int end); glm::vec3 classifyFieldline(FieldlineEnd fEnd, FieldlineEnd bEnd); diff --git a/src/rendering/renderablefieldlines.cpp b/src/rendering/renderablefieldlines.cpp index 64dab8daaa..4527ad2cc6 100644 --- a/src/rendering/renderablefieldlines.cpp +++ b/src/rendering/renderablefieldlines.cpp @@ -94,14 +94,6 @@ RenderableFieldlines::RenderableFieldlines(const ghoul::Dictionary& dictionary) } } - _seedpointsProgram = new ghoul::opengl::ProgramObject("SeedpointsProgram"); - ghoul::opengl::ShaderObject* seedpointVertexShader = new ghoul::opengl::ShaderObject(ghoul::opengl::ShaderObject::ShaderTypeVertex, - "/home/hhellteg/openspace/openspace-data/scene/fieldlines/seedPoints.vert"); - ghoul::opengl::ShaderObject* seedpointFragmentShader = new ghoul::opengl::ShaderObject(ghoul::opengl::ShaderObject::ShaderTypeFragment, - "/home/hhellteg/openspace/openspace-data/scene/fieldlines/seedPoints.frag"); - _seedpointsProgram->attachObject(seedpointVertexShader); - _seedpointsProgram->attachObject(seedpointFragmentShader); - if(dictionary.hasKey("UpdateOnSave")) { dictionary.getValue("UpdateOnSave", _programUpdateOnSave); } @@ -129,14 +121,17 @@ bool RenderableFieldlines::initialize() { _lineCount.push_back(fieldlinesData[j].size()/2.0); prevEnd = prevEnd + fieldlinesData[j].size()/2.0; - _seedpointStart.push_back(j); - _seedpointCount.push_back(1); - vertexData.insert( vertexData.end(), fieldlinesData[j].begin(), fieldlinesData[j].end()); } } + LDEBUG("Number of vertices : " << vertexData.size()/2.0); - LDEBUG("Vertex orginal : " << vertexData.size()/2.0); + // Give seedpoints a color for visualizing as GL_POINTS + std::vector seedPointsData; + for (int i = 0; i < _seedPoints.size(); ++i) { + seedPointsData.push_back(_seedPoints[i]); + seedPointsData.push_back(glm::vec3(1.0, 0.0, 1.0)); + } // ------ FIELDLINES ----------------- GLuint vertexPositionBuffer; @@ -159,22 +154,27 @@ bool RenderableFieldlines::initialize() { glBindBuffer(GL_ARRAY_BUFFER, 0); //unbind buffer glBindVertexArray(0); //unbind array -// // ------ SEEDPOINTS ----------------- + // ------ SEEDPOINTS ----------------- GLuint seedpointPositionBuffer; glGenVertexArrays(1, &_seedpointVAO); // generate array glBindVertexArray(_seedpointVAO); // bind array glGenBuffers(1, &seedpointPositionBuffer); // generate buffer glBindBuffer(GL_ARRAY_BUFFER, seedpointPositionBuffer); // bind buffer - glBufferData(GL_ARRAY_BUFFER, _seedPoints.size()*sizeof(glm::vec3), &_seedPoints.front(), GL_STATIC_DRAW); + glBufferData(GL_ARRAY_BUFFER, seedPointsData.size()*sizeof(glm::vec3), &seedPointsData.front(), GL_STATIC_DRAW); // Vertex positions - GLuint seedpointLocation = 0; - glEnableVertexAttribArray(seedpointLocation); - glVertexAttribPointer(seedpointLocation, 3, GL_FLOAT, GL_FALSE, sizeof(glm::vec3), reinterpret_cast(0)); + glEnableVertexAttribArray(vertexLocation); + glVertexAttribPointer(vertexLocation, 3, GL_FLOAT, GL_FALSE, 2*sizeof(glm::vec3), reinterpret_cast(0)); + + // Texture coordinates + glEnableVertexAttribArray(texcoordLocation); + glVertexAttribPointer(texcoordLocation, 3, GL_FLOAT, GL_FALSE, 2*sizeof(glm::vec3), (void*)(sizeof(glm::vec3))); glBindBuffer(GL_ARRAY_BUFFER, 0); //unbind buffer glBindVertexArray(0); //unbind array + glPointSize(5); // size of seedpoints + // ------ SETUP SHADERS ----------------- auto privateCallback = [this](const ghoul::filesystem::File& file) { safeShaderCompilation(); @@ -187,9 +187,6 @@ bool RenderableFieldlines::initialize() { _fieldlinesProgram->compileShaderObjects(); _fieldlinesProgram->linkProgramObject(); - _seedpointsProgram->compileShaderObjects(); - _seedpointsProgram->linkProgramObject(); - return true; } @@ -211,25 +208,15 @@ void RenderableFieldlines::render(const Camera* camera, const psc& thisPosition) _shaderMutex->lock(); _fieldlinesProgram->activate(); _fieldlinesProgram->setUniform("modelViewProjection", transform); - glBindVertexArray(_VAO); glMultiDrawArrays(GL_LINE_STRIP, &_lineStart[0], &_lineCount[0], _lineStart.size()); - glBindVertexArray(0); - - _fieldlinesProgram->deactivate(); - _shaderMutex->unlock(); // ------ SEEDPOINTS ----------------- - _shaderMutex->lock(); - _seedpointsProgram->activate(); - _seedpointsProgram->setUniform("modelViewProjection", transform); - glBindVertexArray(_seedpointVAO); - glPointSize(5); glMultiDrawArrays(GL_POINTS, &_lineStart[0], &_lineCount[0], _seedPoints.size()); glBindVertexArray(0); - _seedpointsProgram->deactivate(); + _fieldlinesProgram->deactivate(); _shaderMutex->unlock(); } @@ -251,7 +238,11 @@ std::vector > RenderableFieldlines::getFieldlinesData(std KameleonWrapper::Model model; std::vector > fieldlinesData; + bool classification = false, lorentz = false; + glm::vec3 fieldlineColor = glm::vec3(1.0, 1.0, 1.0); // default color if no color or classification is specified + if (hintsDictionary.hasKey("Model") && hintsDictionary.getValue("Model", modelString)) { + // ------ MODEL ----------------- if (modelString == "BATSRUS") { model = KameleonWrapper::Model::BATSRUS; } else if (modelString == "ENLIL") { @@ -262,25 +253,35 @@ std::vector > RenderableFieldlines::getFieldlinesData(std return fieldlinesData; } + // ------ VARIBLES / LORENTZ ----------------- if (hintsDictionary.hasKey("Variables")) { bool xVar, yVar, zVar; xVar = hintsDictionary.getValue("Variables.1", xVariable); - yVar = hintsDictionary.getValue("Variables.2", yVariable); - zVar = hintsDictionary.getValue("Variables.3", zVariable); + if (xVar && xVariable == "Lorentz") { + lorentz = true; + } else { - if (!xVar || !yVar || !zVar) { - LWARNING("Error reading variables! Must be 3 and must exist in CDF data"); - return fieldlinesData; + yVar = hintsDictionary.getValue("Variables.2", yVariable); + zVar = hintsDictionary.getValue("Variables.3", zVariable); + + if (!xVar || !yVar || !zVar) { + LWARNING("Error reading variables! Must be 3 and must exist in CDF data"); + return fieldlinesData; + } } } else { LWARNING("Hints does not specify valid 'Variables'"); return fieldlinesData; } + // ------ STEPSIZE ----------------- if (!hintsDictionary.hasKey("Stepsize") || !hintsDictionary.getValue("Stepsize", stepSize)) { LDEBUG("No stepsize set for fieldlines. Setting to default value (" << stepSize << ")"); } + + // ------ SEEDPOINTS --------------- ghoul::Dictionary seedpointsDictionary; + _seedPoints.clear(); if (hintsDictionary.hasKey("Seedpoints") && hintsDictionary.getValue("Seedpoints", seedpointsDictionary)) { glm::vec3 seedPos; for (auto index : seedpointsDictionary.keys()) { @@ -288,10 +289,22 @@ std::vector > RenderableFieldlines::getFieldlinesData(std _seedPoints.push_back(seedPos); } } + + // ------ CLASSIFICATION & COLOR ----------- + hintsDictionary.getValue("Color", fieldlineColor); + hintsDictionary.getValue("Classification", classification); } KameleonWrapper kw(filename, model); - fieldlinesData = kw.getFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize); + if (lorentz) { + fieldlinesData = kw.getLorentzTrajectories(_seedPoints, fieldlineColor, stepSize); + } else { + if (classification) + fieldlinesData = kw.getClassifiedFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize); + else + fieldlinesData = kw.getFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize, fieldlineColor); + } + return fieldlinesData; } diff --git a/src/util/kameleonwrapper.cpp b/src/util/kameleonwrapper.cpp index 314074c0a0..ead509d8e0 100644 --- a/src/util/kameleonwrapper.cpp +++ b/src/util/kameleonwrapper.cpp @@ -229,7 +229,7 @@ float* KameleonWrapper::getUniformSampledVectorValues(const std::string& xVar, c return data; } -std::vector > KameleonWrapper::getFieldLines( +std::vector > KameleonWrapper::getClassifiedFieldLines( const std::string& xVar, const std::string& yVar, const std::string& zVar, std::vector seedPoints, float stepSize ) { @@ -268,6 +268,67 @@ std::vector > KameleonWrapper::getFieldLines( return fieldLines; } +std::vector > KameleonWrapper::getFieldLines( + const std::string& xVar, const std::string& yVar, + const std::string& zVar, std::vector seedPoints, + float stepSize, glm::vec3 color ) { + assert(_model && _interpolator); + assert(_type == Model::ENLIL || _type == Model::BATSRUS); + LINFO("Creating " << seedPoints.size() << " fieldlines from variables " << xVar << " " << yVar << " " << zVar); + + std::vector fLine, bLine; + std::vector > fieldLines; + FieldlineEnd forwardEnd, backEnd; + + if (_type == Model::BATSRUS) { + for (glm::vec3 seedPoint : seedPoints) { + fLine = traceCartesianFieldline(xVar, yVar, zVar, seedPoint, stepSize, TraceDirection::FORWARD, forwardEnd); + bLine = traceCartesianFieldline(xVar, yVar, zVar, seedPoint, stepSize, TraceDirection::BACK, backEnd); + + bLine.insert(bLine.begin(), fLine.rbegin(), fLine.rend()); + + // write colors + std::vector line; + for (glm::vec3 position : bLine) { + line.push_back(position); + line.push_back(color); + } + + fieldLines.push_back(line); + } + } else { + LERROR("Fieldlines are only supported for BATSRUS model"); + } + + return fieldLines; +} + +std::vector > KameleonWrapper::getLorentzTrajectories( + std::vector seedPoints, glm::vec3 color, float stepsize) { + LINFO("Creating " << seedPoints.size() << " Lorentz force trajectories"); + + std::vector > trajectories; + std::vector plusTraj, minusTraj; + + for (auto seedPoint : seedPoints) { + plusTraj = traceLorentzTrajectory(seedPoint, stepsize, 1.0); + minusTraj = traceLorentzTrajectory(seedPoint, stepsize, -1.0); + + minusTraj.insert(minusTraj.begin(), plusTraj.rbegin(), plusTraj.rend()); + + // write colors + std::vector trajectory; + for (glm::vec3 position : minusTraj) { + trajectory.push_back(position); + trajectory.push_back(color); + } + trajectories.push_back(trajectory); + + } + + return trajectories; +} + std::vector KameleonWrapper::traceCartesianFieldline( const std::string& xVar, const std::string& yVar, const std::string& zVar, glm::vec3 seedPoint, @@ -276,8 +337,7 @@ std::vector KameleonWrapper::traceCartesianFieldline( glm::vec3 color, pos, k1, k2, k3, k4; std::vector line; float stepX, stepY, stepZ; - int numSteps = 0; - int maxSteps = 5000; + int numSteps = 0, maxSteps = 5000; pos = seedPoint; _model->loadVariable(xVar); @@ -300,7 +360,7 @@ std::vector KameleonWrapper::traceCartesianFieldline( k1.y = _interpolator->interpolate(yID, pos.x, pos.y, pos.z); k1.z = _interpolator->interpolate(zID, pos.x, pos.y, pos.z); k1 = (float)direction*glm::normalize(k1); -// stepX*=stepSize, stepY*=stepSize, stepZ*=stepSize; + stepX=stepX*stepSize, stepY=stepY*stepSize, stepZ=stepZ*stepSize; k2.x = _interpolator->interpolate(xID, pos.x+(stepX/2.0)*k1.x, pos.y+(stepY/2.0)*k1.y, pos.z+(stepZ/2.0)*k1.z); k2.y = _interpolator->interpolate(yID, pos.x+(stepX/2.0)*k1.x, pos.y+(stepY/2.0)*k1.y, pos.z+(stepZ/2.0)*k1.z); k2.z = _interpolator->interpolate(zID, pos.x+(stepX/2.0)*k1.x, pos.y+(stepY/2.0)*k1.y, pos.z+(stepZ/2.0)*k1.z); @@ -335,6 +395,96 @@ std::vector KameleonWrapper::traceCartesianFieldline( return line; } +std::vector KameleonWrapper::traceLorentzTrajectory(glm::vec3 seedPoint, + float stepsize, float eCharge) { + glm::vec3 B, E, v0, k1, k2, k3, k4, sPos, tmpV; + float stepX = stepsize, stepY = stepsize, stepZ = stepsize; + + long int bxID = _model->getVariableID("bx"); + long int byID = _model->getVariableID("by"); + long int bzID = _model->getVariableID("bz"); + long int jxID = _model->getVariableID("jx"); + long int jyID = _model->getVariableID("jy"); + long int jzID = _model->getVariableID("jz"); + + std::vector trajectory; + glm::vec3 pos = seedPoint; + int numSteps = 0, maxSteps = 5000; + v0.x = _interpolator->interpolate("ux", pos.x, pos.y, pos.z); + v0.y = _interpolator->interpolate("uy", pos.x, pos.y, pos.z); + v0.z = _interpolator->interpolate("uz", pos.x, pos.y, pos.z); + v0 = glm::normalize(v0); + + // While we are inside the models boundries and not inside earth + while ((pos.x < _xMax && pos.x > _xMin && pos.y < _yMax && pos.y > _yMin && + pos.z < _zMax && pos.z > _zMin) && !(pos.x*pos.x + pos.y*pos.y + pos.z*pos.z < 1.0)) { + + // Save position + trajectory.push_back(pos); + + // Calculate new position with Lorentz force quation and Runge-Kutta 4th order + B.x = _interpolator->interpolate(bxID, pos.x, pos.y, pos.z); + B.y = _interpolator->interpolate(byID, pos.x, pos.y, pos.z); + B.z = _interpolator->interpolate(bzID, pos.x, pos.y, pos.z); + E.x = _interpolator->interpolate(jxID, pos.x, pos.y, pos.z); + E.y = _interpolator->interpolate(jyID, pos.x, pos.y, pos.z); + E.z = _interpolator->interpolate(jzID, pos.x, pos.y, pos.z); + k1 = eCharge*(E + glm::cross(v0, B)); + k1 = glm::normalize(k1); + + sPos = glm::vec3( pos.x+(stepX/2.0)*v0.x+(stepX*stepX/8.0)*k1.x, + pos.y+(stepY/2.0)*v0.y+(stepY*stepY/8.0)*k1.y, + pos.z+(stepZ/2.0)*v0.z+(stepZ*stepZ/8.0)*k1.z); + B.x = _interpolator->interpolate(bxID, sPos.x, sPos.y, sPos.z); + B.y = _interpolator->interpolate(byID, sPos.x, sPos.y, sPos.z); + B.z = _interpolator->interpolate(bzID, sPos.x, sPos.y, sPos.z); + E.x = _interpolator->interpolate(jxID, sPos.x, sPos.y, sPos.z); + E.y = _interpolator->interpolate(jyID, sPos.x, sPos.y, sPos.z); + E.z = _interpolator->interpolate(jzID, sPos.x, sPos.y, sPos.z); + tmpV = v0+(stepX/2.0f)*k1; + k2 = eCharge*(E + glm::cross(tmpV, B)); + k2 = glm::normalize(k2); + + B.x = _interpolator->interpolate(bxID, sPos.x, sPos.y, sPos.z); + B.y = _interpolator->interpolate(byID, sPos.x, sPos.y, sPos.z); + B.z = _interpolator->interpolate(bzID, sPos.x, sPos.y, sPos.z); + E.x = _interpolator->interpolate(jxID, sPos.x, sPos.y, sPos.z); + E.y = _interpolator->interpolate(jyID, sPos.x, sPos.y, sPos.z); + E.z = _interpolator->interpolate(jzID, sPos.x, sPos.y, sPos.z); + tmpV = v0+(stepX/2.0f)*k2; + k3 = eCharge*(E + glm::cross(tmpV, B)); + k3 = glm::normalize(k3); + + sPos = glm::vec3( pos.x+stepX*v0.x+(stepX*stepX/2.0)*k1.x, + pos.y+stepY*v0.y+(stepY*stepY/2.0)*k1.y, + pos.z+stepZ*v0.z+(stepZ*stepZ/2.0)*k1.z); + B.x = _interpolator->interpolate(bxID, sPos.x, sPos.y, sPos.z); + B.y = _interpolator->interpolate(byID, sPos.x, sPos.y, sPos.z); + B.z = _interpolator->interpolate(bzID, sPos.x, sPos.y, sPos.z); + E.x = _interpolator->interpolate(jxID, sPos.x, sPos.y, sPos.z); + E.y = _interpolator->interpolate(jyID, sPos.x, sPos.y, sPos.z); + E.z = _interpolator->interpolate(jzID, sPos.x, sPos.y, sPos.z); + tmpV = v0+stepX*k3; + k4 = eCharge*(E + glm::cross(tmpV, B)); + k4 = glm::normalize(k4); + + pos.x = pos.x + stepX*v0.x + (stepX*stepX/6.0)*(k1.x + k2.x + k3.x); + pos.y = pos.y + stepY*v0.y + (stepY*stepY/6.0)*(k1.y + k2.y + k3.y); + pos.z = pos.z + stepZ*v0.z + (stepZ*stepZ/6.0)*(k1.z + k2.z + k3.z); + + v0.x = v0.x + (stepX/6.0)*(k1.x + 2.0*k2.x + 2.0*k3.x + k4.z); + v0.y = v0.y + (stepY/6.0)*(k1.y + 2.0*k2.y + 2.0*k3.y + k4.y); + v0.z = v0.z + (stepZ/6.0)*(k1.z + 2.0*k2.z + 2.0*k3.z + k4.z); + + ++numSteps; + if (numSteps > maxSteps) { + LDEBUG("Max number of steps taken (" << maxSteps <<")"); + break; + } + } + return trajectory; +} + void KameleonWrapper::getGridVariables(std::string& x, std::string& y, std::string& z) { // get the grid system string std::string gridSystem = _model->getGlobalAttribute("grid_system_1").getAttributeString();