Skip to content

Commit

Permalink
Use flying edges for mesh generation (#1741)
Browse files Browse the repository at this point in the history 
* adding-new-tables-for-flying-edges

* adding codes for flying-edges algo, still fixing left

* fixing laplace smoothing, only cleanups left

Signed-off-by: Perminder Singh <[email protected]>

---------

Signed-off-by: Perminder <[email protected]>
Signed-off-by: Geoff Hutchison <[email protected]>
Signed-off-by: Perminder Singh <[email protected]>
  • Loading branch information
@perminder-17
perminder-17 authored Nov 16, 2024
1 parent 990956e commit 141745d
Show file tree
Hide file tree
Showing 9 changed files with 1,514 additions and 640 deletions.
156 changes: 153 additions & 3 deletions avogadro/core/cube.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@ bool Cube::setLimits(const Vector3& min_, const Vector3& max_,
m_data.resize(m_points.x() * m_points.y() * m_points.z());
return true;
}

bool Cube::setLimits(const Vector3& min_, const Vector3& max_, float spacing_)
{
Vector3 delta = max_ - min_;
Expand Down Expand Up @@ -194,8 +193,159 @@ float Cube::value(int i, int j, int k) const
return 0.0;
}

float Cube::value(const Vector3i& pos) const
std::array<float, 3> Cube::computeGradient(int i, int j, int k) const
{
int nx = m_points.x();
int ny = m_points.y();
int nz = m_points.z();
int dataIdx = (i * ny * nz) + (j * nz) + k;

std::array<std::array<float, 2>, 3> x;
std::array<float, 3> run;

// X-direction
if (i == 0) {
x[0][0] = m_data[dataIdx + ny * nz];
x[0][1] = m_data[dataIdx];
run[0] = m_spacing.x();
} else if (i == (nx - 1)) {
x[0][0] = m_data[dataIdx];
x[0][1] = m_data[dataIdx - ny * nz];
run[0] = m_spacing.x();
} else {
x[0][0] = m_data[dataIdx + ny * nz];
x[0][1] = m_data[dataIdx - ny * nz];
run[0] = 2 * m_spacing.x();
}

// Y-direction
if (j == 0) {
x[1][0] = m_data[dataIdx + nz];
x[1][1] = m_data[dataIdx];
run[1] = m_spacing.y();
} else if (j == (ny - 1)) {
x[1][0] = m_data[dataIdx];
x[1][1] = m_data[dataIdx - nz];
run[1] = m_spacing.y();
} else {
x[1][0] = m_data[dataIdx + nz];
x[1][1] = m_data[dataIdx - nz];
run[1] = 2 * m_spacing.y();
}

// Z-direction
if (k == 0) {
x[2][0] = m_data[dataIdx + 1];
x[2][1] = m_data[dataIdx];
run[2] = m_spacing.z();
} else if (k == (nz - 1)) {
x[2][0] = m_data[dataIdx];
x[2][1] = m_data[dataIdx - 1];
run[2] = m_spacing.z();
} else {
x[2][0] = m_data[dataIdx + 1];
x[2][1] = m_data[dataIdx - 1];
run[2] = 2 * m_spacing.z();
}

std::array<float, 3> ret;

ret[0] = (x[0][1] - x[0][0]) / run[0];
ret[1] = (x[1][1] - x[1][0]) / run[1];
ret[2] = (x[2][1] - x[2][0]) / run[2];

return ret;
}

std::array<std::array<float, 3>, 8>
Cube::getGradCube(int i, int j, int k) const
{
std::array<std::array<float, 3>, 8> grad;

grad[0] = computeGradient(i, j, k);
grad[1] = computeGradient(i + 1, j, k);
grad[2] = computeGradient(i + 1, j + 1, k);
grad[3] = computeGradient(i, j + 1, k);
grad[4] = computeGradient(i, j, k + 1);
grad[5] = computeGradient(i + 1, j, k + 1);
grad[6] = computeGradient(i + 1, j + 1, k + 1);
grad[7] = computeGradient(i, j + 1, k + 1);

return grad;
}

std::array<float, 8> Cube::getValsCube(int i, int j, int k) const
{
std::array<float, 8> vals;

vals[0] = getData(i, j, k);
vals[1] = getData(i + 1, j, k);
vals[2] = getData(i + 1, j + 1, k);
vals[3] = getData(i, j + 1, k);
vals[4] = getData(i, j, k + 1);
vals[5] = getData(i + 1, j, k + 1);
vals[6] = getData(i + 1, j + 1, k + 1);
vals[7] = getData(i, j + 1, k + 1);

return vals;
}


float Cube::getData(int i, int j, int k) const
{
int nx = m_points.x();
int ny = m_points.y();
int nz = m_points.z();
return m_data[(i * ny * nz) + (j * nz) + k];
}


std::array<std::array<float, 3>, 8> Cube::getPosCube(int i, int j, int k) const
{

std::array<std::array<float, 3>, 8> pos;

float xpos = m_min.x() + (i * m_spacing.x());
float ypos = m_min.y() + (j * m_spacing.y());
float zpos = m_min.z() + (k * m_spacing.z());

pos[0][0] = xpos;
pos[0][1] = ypos;
pos[0][2] = zpos;

pos[1][0] = xpos + m_spacing.x();
pos[1][1] = ypos;
pos[1][2] = zpos;

pos[2][0] = xpos + m_spacing.x();
pos[2][1] = ypos + m_spacing.y();
pos[2][2] = zpos;

pos[3][0] = xpos;
pos[3][1] = ypos + m_spacing.y();
pos[3][2] = zpos;

pos[4][0] = xpos;
pos[4][1] = ypos;
pos[4][2] = zpos + m_spacing.z();

pos[5][0] = xpos + m_spacing.x();
pos[5][1] = ypos;
pos[5][2] = zpos + m_spacing.z();

pos[6][0] = xpos + m_spacing.x();
pos[6][1] = ypos + m_spacing.y();
pos[6][2] = zpos + m_spacing.z();

pos[7][0] = xpos;
pos[7][1] = ypos + m_spacing.y();
pos[7][2] = zpos + m_spacing.z();

return pos;
}

float Cube::value(const Vector3i& pos) const
{
unsigned int index =
pos.x() * m_points.y() * m_points.z() + pos.y() * m_points.z() + pos.z();
if (index < m_data.size())
Expand Down Expand Up @@ -294,4 +444,4 @@ bool Cube::fillStripe(
return true;
}

} // End Avogadro namespace
} // End Avogadro namespace
65 changes: 61 additions & 4 deletions avogadro/core/cube.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "vector.h"

#include <vector>
#include <array>

namespace Avogadro::Core {

Expand All @@ -21,6 +22,7 @@ class Mutex;
* @class Cube cube.h <avogadro/core/cube.h>
* @brief Provide a data structure for regularly spaced 3D grids.
* @author Marcus D. Hanwell
* @author Perminder Singh
*/

class AVOGADROCORE_EXPORT Cube
Expand Down Expand Up @@ -115,7 +117,7 @@ class AVOGADROCORE_EXPORT Cube
bool setLimits(const Molecule& mol, float spacing, float padding);

/**
* @return Vector containing all the data in a one-dimensional array.
* @return Pointer to the vector containing all the data in a one-dimensional array.
*/
std::vector<float>* data();
const std::vector<float>* data() const;
Expand All @@ -130,6 +132,11 @@ class AVOGADROCORE_EXPORT Cube
*/
bool addData(const std::vector<float>& values);

/**
* @return Const iterator to the beginning of a specific row in the cube data.
*/
std::vector<float>::const_iterator getRowIter(int j, int k) const;

/**
* @return Index of the point closest to the position supplied.
* @param pos Position to get closest index for.
Expand Down Expand Up @@ -198,7 +205,7 @@ class AVOGADROCORE_EXPORT Cube
* @param value Value to fill the cube with.
*/
void fill(float value);

/**
* Sets all indices in a Z stripe of the cube to the specified value.
* @param i x component of the position.
Expand All @@ -212,12 +219,12 @@ class AVOGADROCORE_EXPORT Cube
);

/**
* @return The minimum value at any point in the Cube.
* @return The minimum value at any point in the Cube.
*/
float minValue() const { return m_minValue; }

/**
* @return The maximum value at any point in the Cube.
* @return The maximum value at any point in the Cube.
*/
float maxValue() const { return m_maxValue; }

Expand All @@ -232,6 +239,56 @@ class AVOGADROCORE_EXPORT Cube
*/
Mutex* lock() const { return m_lock; }

/**
* Compute the gradient at a specific point in the cube.
* @param i x index
* @param j y index
* @param k z index
* @return Gradient vector at the specified point.
*/
std::array<float, 3> computeGradient(int i, int j, int k) const;

/**
* Get the values of the eight corners of a cube defined by the indices (i, j, k).
* @param i x index
* @param j y index
* @param k z index
* @return Array of values at the eight corners.
*/
std::array<float, 8> getValsCube(int i, int j, int k) const;

/**
* Get the gradients at the eight corners of the cube defined by the indices (i, j, k).
* @param i x index
* @param j y index
* @param k z index
* @return Array of gradients at the eight corners. Each gradient is a 3D vector.
*/
std::array<std::array<float, 3>, 8> getGradCube(int i, int j, int k) const;

/**
* Get the data value at the specified indices.
* @param i x index
* @param j y index
* @param k z index
* @return Value at the specified indices.
*/
float getData(int i, int j, int k) const;

/**
* Retrieves the positions of the eight corners of a cube at grid indices (i, j, k).
*
* The indices (i, j, k) are converted to real-space positions (xpos, ypos, zpos), mapping grid indices to physical coordinates.
* The method returns a cube that spans one step in each of the x, y, and z directions, with step sizes defined by `m_spacing`.
*
* @param i X-index.
* @param j Y-index.
* @param k Z-index.
* @return A `std::array` of eight `(x, y, z)` coordinates representing the cube's corners.
*/

std::array<std::array<float, 3>, 8> getPosCube(int i, int j, int k) const;

protected:
std::vector<float> m_data;
Vector3 m_min, m_max, m_spacing;
Expand Down
Loading

0 comments on commit 141745d

Please sign in to comment.