Added support for multiple fieldlines in a single module. Stepsize is now adaptive based on the local grid size in the model

This commit is contained in:
Hans-Christian Helltegen
2014-06-16 11:39:10 -04:00
parent 6f40de5f37
commit 56c5c259fb
4 changed files with 132 additions and 125 deletions

View File

@@ -52,8 +52,10 @@ public:
virtual void update();
private:
ghoul::Dictionary _hintsDictionary;
std::string _filename;
std::vector<std::vector<glm::vec3> > getFieldlinesData(std::string filename, ghoul::Dictionary hintsDictionary);
std::vector<ghoul::Dictionary> _hintsDictionaries;
std::vector<std::string> _filenames;
std::vector<glm::vec3> _seedPoints;
ghoul::opengl::ProgramObject *_fieldlinesProgram, *_seedpointsProgram;
@@ -69,8 +71,6 @@ private:
bool _programUpdateOnSave;
void safeShaderCompilation();
};
} // namespace openspace

View File

@@ -49,23 +49,15 @@ RenderableFieldlines::RenderableFieldlines(const ghoul::Dictionary& dictionary)
file = findPath(file);
if (file != "") {
// parse hints
// read hints into dictionary
ghoul::Dictionary hintsDictionary;
if(fieldline.hasKey("Hints"))
fieldline.getValue("Hints", hintsDictionary);
// TODO Vectors of filenames and dictionaries
_filename = file;
_hintsDictionary = hintsDictionary;
_filenames.push_back(file);
_hintsDictionaries.push_back(hintsDictionary);
ghoul::Dictionary seedpointsDictionary;
if (fieldline.hasKey("Seedpoints") && fieldline.getValue("Seedpoints", seedpointsDictionary)) {
glm::vec3 tmpVal;
for (int i = 0; i < seedpointsDictionary.keys().size(); ++i) {
fieldline.getValue("Seedpoints."+std::to_string(i+1), tmpVal);
_seedPoints.push_back(tmpVal);
}
}
} else
LERROR("File not found!");
}
@@ -122,65 +114,30 @@ RenderableFieldlines::~RenderableFieldlines() {
}
bool RenderableFieldlines::initialize() {
assert(_filename != "");
assert(_hintsDictionary.size() != 0);
assert(_seedPoints.size() != 0);
assert(_filenames.size() != 0);
assert(_hintsDictionaries.size() != 0);
std::string modelString;
std::vector<std::vector<glm::vec3> > fieldlinesData;
float stepSize;
std::string xVariable, yVariable, zVariable;
KameleonWrapper::Model model;
if (_hintsDictionary.hasKey("Model") && _hintsDictionary.getValue("Model", modelString)) {
if (modelString == "BATSRUS") {
model = KameleonWrapper::Model::BATSRUS;
} else if (modelString == "ENLIL") {
LWARNING("ENLIL model not supported for fieldlines");
return false;
} else {
LWARNING("Hints does not specify a valid 'Model'");
return false;
}
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 || !yVar || !zVar) {
LWARNING("Error reading variables! Must be 3 and must exist in CDF data");
return false;
}
} else {
LWARNING("Hints does not specify valid 'Variables'");
return false;
}
if (!_hintsDictionary.hasKey("Stepsize") || !_hintsDictionary.getValue("Stepsize", stepSize)) {
LDEBUG("No stepsize set for fieldlines. Setting to default value (0.5)");
stepSize = 0.5;
}
}
KameleonWrapper kw(_filename, model);
fieldlinesData = kw.getFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize);
std::vector<glm::vec3> vertexData;
int prevEnd = 0;
std::vector<glm::vec3> vertexData;
std::vector<std::vector<glm::vec3> > fieldlinesData;
for (int i = 0; i < fieldlinesData.size(); i++) {
_lineStart.push_back(prevEnd);
_lineCount.push_back(fieldlinesData[i].size()/2.0);
prevEnd = prevEnd + fieldlinesData[i].size()/2.0;
for (int i = 0; i < _filenames.size(); ++i) {
fieldlinesData = getFieldlinesData(_filenames[i], _hintsDictionaries[i]);
_seedpointStart.push_back(i);
_seedpointCount.push_back(1);
for (int j = 0; j < fieldlinesData.size(); j++) {
_lineStart.push_back(prevEnd);
_lineCount.push_back(fieldlinesData[j].size()/2.0);
prevEnd = prevEnd + fieldlinesData[j].size()/2.0;
vertexData.insert( vertexData.end(), fieldlinesData[i].begin(), fieldlinesData[i].end());
_seedpointStart.push_back(j);
_seedpointCount.push_back(1);
vertexData.insert( vertexData.end(), fieldlinesData[j].begin(), fieldlinesData[j].end());
}
}
LDEBUG("Vertex orginal : " << vertexData.size()/2.0);
// ------ FIELDLINES -----------------
GLuint vertexPositionBuffer;
glGenVertexArrays(1, &_VAO); // generate array
@@ -203,20 +160,20 @@ bool RenderableFieldlines::initialize() {
glBindVertexArray(0); //unbind array
// // ------ 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);
//
// // Vertex positions
// GLuint seedpointLocation = 0;
// glEnableVertexAttribArray(seedpointLocation);
// glVertexAttribPointer(seedpointLocation, 3, GL_FLOAT, GL_FALSE, 3*sizeof(GLfloat), reinterpret_cast<void*>(0));
//
// glBindBuffer(GL_ARRAY_BUFFER, 0); //unbind buffer
// glBindVertexArray(0); //unbind array
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);
// Vertex positions
GLuint seedpointLocation = 0;
glEnableVertexAttribArray(seedpointLocation);
glVertexAttribPointer(seedpointLocation, 3, GL_FLOAT, GL_FALSE, sizeof(glm::vec3), reinterpret_cast<void*>(0));
glBindBuffer(GL_ARRAY_BUFFER, 0); //unbind buffer
glBindVertexArray(0); //unbind array
// ------ SETUP SHADERS -----------------
auto privateCallback = [this](const ghoul::filesystem::File& file) {
@@ -247,8 +204,8 @@ void RenderableFieldlines::render(const Camera* camera, const psc& thisPosition)
transform = transform*camTransform;
transform = glm::translate(transform, relative.vec3());
transform = glm::rotate(transform, -90.0f, glm::vec3(1.0, 0.0, 0.0));
transform = glm::scale(transform, glm::vec3(0.1));
transform = glm::rotate(transform, -90.0f, glm::vec3(1.0, 0.0, 0.0)); // Model has positive Z as up
transform = glm::scale(transform, glm::vec3(0.1)); // Scale to avoid depth buffer problems
// ------ FIELDLINES -----------------
_shaderMutex->lock();
@@ -262,18 +219,18 @@ void RenderableFieldlines::render(const Camera* camera, const psc& thisPosition)
_fieldlinesProgram->deactivate();
_shaderMutex->unlock();
// // ------ SEEDPOINTS -----------------
// _shaderMutex->lock();
// _seedpointsProgram->activate();
// _seedpointsProgram->setUniform("modelViewProjection", transform);
//
// glBindVertexArray(_seedpointVAO);
// glPointSize(3);
// glMultiDrawArrays(GL_POINTS, &_seedpointStart[0], &_seedpointCount[0], _seedPoints.size());
// glBindVertexArray(0);
//
// _seedpointsProgram->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();
_shaderMutex->unlock();
}
void RenderableFieldlines::update() {
@@ -287,4 +244,55 @@ void RenderableFieldlines::safeShaderCompilation() {
_shaderMutex->unlock();
}
std::vector<std::vector<glm::vec3> > RenderableFieldlines::getFieldlinesData(std::string filename, ghoul::Dictionary hintsDictionary) {
std::string modelString;
float stepSize = 0.5; // default if no stepsize is specified in hints
std::string xVariable, yVariable, zVariable;
KameleonWrapper::Model model;
std::vector<std::vector<glm::vec3> > fieldlinesData;
if (hintsDictionary.hasKey("Model") && hintsDictionary.getValue("Model", modelString)) {
if (modelString == "BATSRUS") {
model = KameleonWrapper::Model::BATSRUS;
} else if (modelString == "ENLIL") {
LWARNING("ENLIL model not supported for fieldlines");
return fieldlinesData;
} else {
LWARNING("Hints does not specify a valid 'Model'");
return fieldlinesData;
}
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 || !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;
}
if (!hintsDictionary.hasKey("Stepsize") || !hintsDictionary.getValue("Stepsize", stepSize)) {
LDEBUG("No stepsize set for fieldlines. Setting to default value (" << stepSize << ")");
}
ghoul::Dictionary seedpointsDictionary;
if (hintsDictionary.hasKey("Seedpoints") && hintsDictionary.getValue("Seedpoints", seedpointsDictionary)) {
glm::vec3 seedPos;
for (auto index : seedpointsDictionary.keys()) {
hintsDictionary.getValue("Seedpoints."+index, seedPos);
_seedPoints.push_back(seedPos);
}
}
}
KameleonWrapper kw(filename, model);
fieldlinesData = kw.getFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize);
return fieldlinesData;
}
} // namespace openspace

View File

@@ -63,8 +63,6 @@ KameleonWrapper::KameleonWrapper(const std::string& filename, Model model): _typ
}
getGridVariables(_xCoordVar, _yCoordVar, _zCoordVar);
LDEBUG("Using coordinate system variables: " << _xCoordVar << ", " << _yCoordVar << ", " << _zCoordVar);
_xMin = _model->getVariableAttribute(_xCoordVar, "actual_min").getAttributeFloat();
_xMax = _model->getVariableAttribute(_xCoordVar, "actual_max").getAttributeFloat();
_yMin = _model->getVariableAttribute(_yCoordVar, "actual_min").getAttributeFloat();
@@ -72,13 +70,6 @@ KameleonWrapper::KameleonWrapper(const std::string& filename, Model model): _typ
_zMin = _model->getVariableAttribute(_zCoordVar, "actual_min").getAttributeFloat();
_zMax = _model->getVariableAttribute(_zCoordVar, "actual_max").getAttributeFloat();
LDEBUG(_xCoordVar << "Min: " << _xMin);
LDEBUG(_xCoordVar << "Max: " << _xMax);
LDEBUG(_yCoordVar << "Min: " << _yMin);
LDEBUG(_yCoordVar << "Max: " << _yMax);
LDEBUG(_zCoordVar << "Min: " << _zMin);
LDEBUG(_zCoordVar << "Max: " << _zMax);
_lastiProgress = -1; // For progressbar
}
@@ -255,6 +246,7 @@ std::vector<std::vector<glm::vec3> > KameleonWrapper::getFieldLines(
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());
// classify
@@ -283,37 +275,43 @@ std::vector<glm::vec3> KameleonWrapper::traceCartesianFieldline(
glm::vec3 color, pos, k1, k2, k3, k4;
std::vector<glm::vec3> line;
float stepX = stepSize, stepY = stepSize, stepZ = stepSize; // Should I do different stepsizes?
float stepX, stepY, stepZ;
int numSteps = 0;
int maxSteps = 5000;
pos = seedPoint;
_model->loadVariable(xVar);
_model->loadVariable(yVar);
_model->loadVariable(zVar);
long int xID = _model->getVariableID(xVar);
long int yID = _model->getVariableID(yVar);
long int zID = _model->getVariableID(zVar);
// 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 < 1.0 && pos.x > -1.0 &&
pos.y < 1.0 && pos.y > -1.0 && pos.z < 1.0 && pos.z > -1.0)) {
pos.z < _zMax && pos.z > _zMin) && !(pos.x*pos.x + pos.y*pos.y + pos.z*pos.z < 1.0)) {
// Save position
line.push_back(pos);
// Calculate new position with Runge-Kutta 4th order
k1.x = _interpolator->interpolate(xVar, pos.x, pos.y, pos.z);
k1.y = _interpolator->interpolate(yVar, pos.x, pos.y, pos.z);
k1.z = _interpolator->interpolate(zVar, pos.x, pos.y, pos.z);
k1.x = _interpolator->interpolate(xID, pos.x, pos.y, pos.z, stepX, stepY, stepZ);
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);
k2.x = _interpolator->interpolate(xVar, 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(yVar, 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(zVar, pos.x+(stepX/2.0)*k1.x, pos.y+(stepY/2.0)*k1.y, pos.z+(stepZ/2.0)*k1.z);
// stepX*=stepSize, stepY*=stepSize, 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);
k2 = (float)direction*glm::normalize(k2);
k3.x = _interpolator->interpolate(xVar, pos.x+(stepX/2.0)*k2.x, pos.y+(stepY/2.0)*k2.y, pos.z+(stepZ/2.0)*k2.z);
k3.y = _interpolator->interpolate(yVar, pos.x+(stepX/2.0)*k2.x, pos.y+(stepY/2.0)*k2.y, pos.z+(stepZ/2.0)*k2.z);
k3.z = _interpolator->interpolate(zVar, pos.x+(stepX/2.0)*k2.x, pos.y+(stepY/2.0)*k2.y, pos.z+(stepZ/2.0)*k2.z);
k3.x = _interpolator->interpolate(xID, pos.x+(stepX/2.0)*k2.x, pos.y+(stepY/2.0)*k2.y, pos.z+(stepZ/2.0)*k2.z);
k3.y = _interpolator->interpolate(yID, pos.x+(stepX/2.0)*k2.x, pos.y+(stepY/2.0)*k2.y, pos.z+(stepZ/2.0)*k2.z);
k3.z = _interpolator->interpolate(zID, pos.x+(stepX/2.0)*k2.x, pos.y+(stepY/2.0)*k2.y, pos.z+(stepZ/2.0)*k2.z);
k3 = (float)direction*glm::normalize(k3);
k4.x = _interpolator->interpolate(xVar, pos.x+stepX*k3.x, pos.y+stepY*k3.y, pos.z+stepZ*k3.z);
k4.y = _interpolator->interpolate(yVar, pos.x+stepX*k3.x, pos.y+stepY*k3.y, pos.z+stepZ*k3.z);
k4.z = _interpolator->interpolate(zVar, pos.x+stepX*k3.x, pos.y+stepY*k3.y, pos.z+stepZ*k3.z);
k4.x = _interpolator->interpolate(xID, pos.x+stepX*k3.x, pos.y+stepY*k3.y, pos.z+stepZ*k3.z);
k4.y = _interpolator->interpolate(yID, pos.x+stepX*k3.x, pos.y+stepY*k3.y, pos.z+stepZ*k3.z);
k4.z = _interpolator->interpolate(zID, pos.x+stepX*k3.x, pos.y+stepY*k3.y, pos.z+stepZ*k3.z);
k4 = (float)direction*glm::normalize(k4);
pos.x = pos.x + (stepX/6.0)*(k1.x + 2.0*k2.x + 2.0*k3.x + k4.x);
pos.y = pos.y + (stepY/6.0)*(k1.y + 2.0*k2.y + 2.0*k3.y + k4.y);
@@ -321,14 +319,15 @@ std::vector<glm::vec3> KameleonWrapper::traceCartesianFieldline(
++numSteps;
if (numSteps > maxSteps) {
LDEBUG("Max number of steps taken");
LDEBUG("Max number of steps taken (" << maxSteps <<")");
break;
}
}
line.push_back(pos);
if (pos.z > 0.0 && pos.x < 1.0 && pos.x > -1.0 && pos.y < 1.0 && pos.y > -1.0)
if (pos.z > 0.0 && (pos.x*pos.x + pos.y*pos.y + pos.z*pos.z < 1.0))
end = FieldlineEnd::NORTH;
else if (pos.z < 0.0 && pos.x < 1.0 && pos.x > -1.0 && pos.y < 1.0 && pos.y > -1.0)
else if (pos.z < 0.0 && (pos.x*pos.x + pos.y*pos.y + pos.z*pos.z < 1.0))
end = FieldlineEnd::SOUTH;
else
end = FieldlineEnd::OUT;