Skip to content

Mesh for anisotropic particles #52

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: refactor-aniso
Choose a base branch
from
Open

Mesh for anisotropic particles #52

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
58fa042
Extend center, lower_bound, upper_bound, shape, step and bin to highe…
mayukh33 Apr 4, 2025
e3dabd9
Fix multi-dimensional indexing
mayukh33 Apr 15, 2025
f374be9
Apply fix and add documentation
mayukh33 Apr 16, 2025
1ece6ab
Minor fix
mayukh33 Apr 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
193 changes: 142 additions & 51 deletions include/flyft/mesh.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,119 +8,210 @@

namespace flyft
{

//! Multidimensional mesh
/*!
* Multidimensional mesh maps N-dimensional indexing for
* the anistropic particles into one-dimensional flux.
*/
class Mesh
{
public:
//! No default constructor
Mesh() = delete;

//! Constructor
/*!
* \param lower_bound Lower bound of the N-dimensional mesh.
* \param upper_bound Upper bound of the N-dimensional mesh.
* \param shape Shape of the index array.
* \param lower_bc Boundary condition of the lower bound.
* \param upper_bc Boundary condition of the upper bound.
*/
Mesh(double lower_bound,
double upper_bound,
int shape,
BoundaryType lower_bc,
BoundaryType upper_bc);

std::shared_ptr<Mesh> slice(int start, int end) const;
//! Slice multidimensional index
/*!
* \param start Multidimensional start index.
* \param end Multidimensional end index.
* \return One-dimensional index.
*
*/
std::shared_ptr<Mesh> slice(const std::vector<int>& start, const std::vector<int>& end) const;

//! Get position on the mesh, defined as center of bin
double center(int i) const;
/*!
* \param i Multidimensional bin index
* \return Center of the multidimensional bin
*/
std::vector<double> center(const std::vector<int>& i) const;

//! Lower bound of entire mesh
double lower_bound() const;
/*!
* \return Lower bound of the multidimensional mesh
*/
std::vector<double> lower_bound() const;

//! Get lower bound of bin
double lower_bound(int i) const;
/*!
* \param i Multidimensional bin index
* \return Lower bound of the multidimensional bin
*/
std::vector<double> lower_bound(const std::vector<int>& i) const;

//! Upper bound of entire mesh
double upper_bound() const;
/*!
* \return Upper bound of the multidimensional mesh
*/
std::vector<double> upper_bound() const;

//! Get upper bound of bin
double upper_bound(int i) const;
/*!
* \param i Multidimensional bin index
* \return Upper bound of the multidimensional bin
*/
std::vector<double> upper_bound(const std::vector<int>& i) const;

//! Get surface area of lower edge of bin
/*!
* \param i Multidimensional bin index
* \return Upper bound of the multidimensional bin
*/
virtual double area(int i) const = 0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For this method, what area is this going to return? Should it return the areas of all lower faces in some defined order? Or, should it take both an index to a bin and an index to a dimension and return the area of that lower face?

Or perhaps, it needs three arguments: bin, face, and +/- direction?


//! Get volume of mesh
/*!
* \return Volume of the mesh
*/
virtual double volume() const = 0;

//! Get volume of bin
/*!
* \param i Multidimensional bin index
* \return Volume of the bin
*/
virtual double volume(int i) const = 0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This needs to be a vector


//! Get the bin for a coordinate
int bin(double x) const;
/*!
* \param x Coordinate position of the bin
* \return Bin index
*/
std::vector<int> bin(double x) const;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is going to need to take a vector, right?


//! Length of the mesh
double L() const;
/*!
* \return Length of the mesh
*/
std::vector<double> L() const;

//! Shape of the mesh
int shape() const;
std::vector<int> shape() const;

//! Step size of the mesh
double step() const;
std::vector<double> step() const;

//! Boundary condition on lower bound of mesh
BoundaryType lower_boundary_condition() const;

//! Boundary condition on upper bound of mesh
BoundaryType upper_boundary_condition() const;

int asShape(double dx) const;

double asLength(int shape) const;

double integrateSurface(int idx, double j_lo, double j_hi) const;
double integrateSurface(int idx, const DataView<double>& j) const;
double integrateSurface(int idx, const DataView<const double>& j) const;

double integrateVolume(int idx, double f) const;
double integrateVolume(int idx, const DataView<double>& f) const;
double integrateVolume(int idx, const DataView<const double>& f) const;

//! Get the start index of the mesh
/*!
* \param dx Step size of the mesh
* \return Start index of the mesh
*/
std::vector<int> asShape(const std::vector<double>& dx) const;

//! Get the length of the mesh
/*!
* \param shape Shape of the mesh
* \return Length of the mesh
*/
std::vector<double> asLength(const std::vector<int>& shape) const;

//! Get the integral of the surface over the mesh
/*!
* \param shape Shape of the mesh
* \return Length of the mesh
*/
double integrateSurface(const std::vector<int>& idx, double j_lo, double j_hi) const;
double integrateSurface(const std::vector<int>& idx, const DataView<double>& j) const;
double integrateSurface(const std::vector<int>& idx, const DataView<const double>& j) const;
Comment on lines +137 to +144
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you think about how these are going to work in the new code? Is it going to integrate over a specific face, or all of them? And what will it integrate?

//! Get the integral of the volume over the mesh
/*!
* \param idx Multidimensional index
* \param f Function to integrate
* \return Integral of the function over the mesh
*/
double integrateVolume(const std::vector<int>& idx, double f) const;
double integrateVolume(const std::vector<int>& idx, const DataView<double>& f) const;
double integrateVolume(const std::vector<int>& idx, const DataView<const double>& f) const;

//! Interpolate a function on the mesh
/*!
* \param x Coordinate position of the bin
* \param f Function to interpolate
* \return Interpolated function value
*
* The function is interpolated using linear interpolation.
*/
template<typename T>
typename std::remove_const<T>::type interpolate(double x, const DataView<T>& f) const;

virtual double gradient(int idx, double f_lo, double f_hi) const = 0;
double gradient(int idx, const DataView<const double>& f) const;
double gradient(int idx, const DataView<double>& f) const;
//! Get the gradient of the mesh
virtual double gradient(const std::vector<int>& idx, double f_lo, double f_hi) const = 0;
double gradient(const std::vector<int>& idx, const DataView<const double>& f) const;
double gradient(const std::vector<int>& idx, const DataView<double>& f) const;

bool operator==(const Mesh& other) const;
bool operator!=(const Mesh& other) const;

protected:
double lower_;
int shape_; //!< Shape of the mesh
BoundaryType lower_bc_;
BoundaryType upper_bc_;
double step_; //!< Spacing between mesh points
int start_;
std::vector<double> lower_;
std::vector<int> shape_; //!< Shape of the mesh
BoundaryType lower_bc_; //!< Boundary condition of the lower bound
BoundaryType upper_bc_; //!< Boundary condition of the upper bound
std::vector<double> step_; //!< Spacing between mesh points
std::vector<int> start_;

void validateBoundaryCondition() const;

virtual std::shared_ptr<Mesh> clone() const = 0;
};

template<typename T>
typename std::remove_const<T>::type Mesh::interpolate(double x, const DataView<T>& f) const
typename std::remove_const<T>::type Mesh::interpolate(std::vector<double> x,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this implementation right? Don't we need to do multidimensional interpolation?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Based on our discussion, this interpolation method probably needs to be updated so that it is operating only on one coordinate, with the rest fixed to a specific bin.

const DataView<T>& f) const
{
const auto idx = bin(x);
const auto x_c = center(idx);
double x_0, x_1;
typename std::remove_const<T>::type f_0, f_1;
if (x < x_c)
{
x_0 = center(idx - 1);
x_1 = x_c;
f_0 = f(idx - 1);
f_1 = f(idx);
}
else
std::vector<T> temp;
for (int k = 0; k < shape_.size(); ++k)
{
x_0 = x_c;
x_1 = center(idx + 1);
f_0 = f(idx);
f_1 = f(idx + 1);
const auto idx = bin(x[k]);
const auto x_c = center(idx);
double x_0, x_1;
typename std::remove_const<T>::type f_0, f_1;
if (x < x_c)
{
x_0 = center(idx - 1);
x_1 = x_c;
f_0 = f(idx - 1);
f_1 = f(idx);
}
else
{
x_0 = x_c;
x_1 = center(idx + 1);
f_0 = f(idx);
f_1 = f(idx + 1);
}
temp.push_back(f_0 + (x - x_0) * (f_1 - f_0) / (x_1 - x_0));
}

return f_0 + (x - x_0) * (f_1 - f_0) / (x_1 - x_0);
return temp;
}

} // namespace flyft
Expand Down
Loading