Skip to content

[Shashwat] Submitting assignment 2 #33

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
shash-sharma wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[Shashwat] Submitting assignment 2 #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
58 changes: 57 additions & 1 deletion src/fd_grad.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,7 @@
#include "fd_grad.h"
#include "fd_partial_derivative.h"
#include <Eigen/Sparse>
#include <iostream>

void fd_grad(
const int nx,
Expand All @@ -8,6 +11,59 @@ void fd_grad(
Eigen::SparseMatrix<double> & G)
{
////////////////////////////////////////////////////////////////////////////
// Add your code here

std::cout << "Building gradient matrix components..." << std::flush;

int dir;

// x-component
dir = 1;
Eigen::SparseMatrix<double> Dx;
fd_partial_derivative(nx, ny, nz, h, dir, Dx);

// y-component
dir = 2;
Eigen::SparseMatrix<double> Dy;
fd_partial_derivative(nx, ny, nz, h, dir, Dy);

// z-component
dir = 3;
Eigen::SparseMatrix<double> Dz;
fd_partial_derivative(nx, ny, nz, h, dir, Dz);

std::cout << "done." << std::endl;

// Concatenate to create gradient
// Note: concatenation method sourced from https://stackoverflow.com/questions/41756428/concatenate-sparse-matrix-eigen

std::cout << "Assembling full gradient matrix..." << std::flush;

G.resize((nx-1)*ny*nz + nx*(ny-1)*nz + nx*ny*(nz-1), nx*ny*nz);

// Using the triplet insertion method based on Eigen's documentation
typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;

for (int ii = 0; ii < G.cols(); ii++)
{
for (Eigen::SparseMatrix<double>::InnerIterator itDx(Dx, ii); itDx; ++itDx)
tripletList.push_back(T(itDx.row(), ii, itDx.value()/h));
// G.insertBack(itDx.row(), ii) = itDx.value()/h;

for (Eigen::SparseMatrix<double>::InnerIterator itDy(Dy, ii); itDy; ++itDy)
tripletList.push_back(T(itDy.row(), ii, itDy.value()/h));
// G.insertBack(itDy.row(), ii) = itDy.value()/h;

for (Eigen::SparseMatrix<double>::InnerIterator itDz(Dz, ii); itDz; ++itDz)
tripletList.push_back(T(itDz.row(), ii, itDz.value()/h));
// G.insertBack(itDz.row(), ii) = itDz.value()/h;
}

G.setFromTriplets(tripletList.begin(), tripletList.end());
G.finalize();

std::cout << "done." << std::endl;


////////////////////////////////////////////////////////////////////////////
}
132 changes: 131 additions & 1 deletion src/fd_interpolate.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
#include "fd_interpolate.h"
#include <iostream>
#include <cmath>

void fd_interpolate(
const int nx,
Expand All @@ -10,6 +12,134 @@ void fd_interpolate(
Eigen::SparseMatrix<double> & W)
{
////////////////////////////////////////////////////////////////////////////
// Add your code here

std::cout << "Building interpolation matrix..." << std::flush;

// ****** DEBUGGING ******
// std::cout << "Rows: " << P.rows() << ", Cols: " << P.cols() << std::endl;
// std::cout << P(0,0) << ", " << P(0,1) << ", " << P(0,2) << std::endl;
// std::cout << corner(0) << ", " << corner(1) << ", " << corner(2) << std::endl;
// ***********************

// Compute total number of grid points
int ng = nx*ny*nz;

// Extract total number of given points
int np = P.rows();

// Since the interpolation matrix multiplies into the list of grid points,
// and interpolates them on to the given points, we know that the size of the
// interpolation matrix must be np x ng. Also, since there are three components,
// we extend the W matrix to concatenate each component on top of each other.

W.resize(3*np, ng);
W.reserve(3*np*4);

// ****** DEBUGGING ******
// std::cout << W.nonZeros() << std::endl;
// ***********************

// For each point in P, we need to find the nearest grid points, and then
// use them to construct W.

// Using the triplet insertion method based on Eigen's documentation
typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;

for (int ii = 0; ii < np; ii++)
{
// Compute distances from this point to the corner
double dist_p_to_corner_x = std::abs(P(ii, 0) - corner(0));
double dist_p_to_corner_y = std::abs(P(ii, 1) - corner(1));
double dist_p_to_corner_z = std::abs(P(ii, 2) - corner(2));

// Divide this distance by grid spacing - this is the number of cells between the point and the corner
double cells_x = dist_p_to_corner_x/h;
double cells_y = dist_p_to_corner_y/h;
double cells_z = dist_p_to_corner_z/h;

// By rounding these numbers up and down, we can figure out the nearest vertices
double x_low = std::floor(cells_x);//*h + corner(0);
double x_high = std::ceil(cells_x);//*h + corner(0);
double y_low = std::floor(cells_y);//*h + corner(1);
double y_high = std::ceil(cells_y);//*h + corner(1);
double z_low = std::floor(cells_z);//*h + corner(2);
double z_high = std::ceil(cells_z);//*h + corner(2);

// Also figure out the *global* indices of the eight surrounding vertices
int i000 = x_low + y_low*nx + z_low*nx*ny;
int i100 = x_high + y_low*nx + z_low*nx*ny;
int i010 = x_low + y_high*nx + z_low*nx*ny;
int i110 = x_high + y_high*nx + z_low*nx*ny;
int i001 = x_low + y_low*nx + z_high*nx*ny;
int i101 = x_high + y_low*nx + z_high*nx*ny;
int i011 = x_low + y_high*nx + z_high*nx*ny;
int i111 = x_high + y_high*nx + z_high*nx*ny;

// ****** DEBUGGING ******
// std::cout << i000 << ", " << i100 << ", " << i010 << ", " << i110 << ", " << i001 << ", " << i101 << ", " << i011 << ", " << i111 << std::endl;
// ***********************

// Now insert weights into appropriate spots of the interpolation matrix
double w_x_low = 1 - (P(ii, 0) - x_low*h + corner(0))/(x_high*h + corner(0) - x_low*h + corner(0));
double w_x_high = (P(ii, 0) - x_low*h + corner(0))/(x_high*h + corner(0) - x_low*h + corner(0));
double w_y_low = 1 - (P(ii, 1) - y_low*h + corner(1))/(y_high*h + corner(1) - y_low*h + corner(1));
double w_y_high = (P(ii, 1) - y_low*h + corner(1))/(y_high*h + corner(1) - y_low*h + corner(1));
double w_z_low = 1 - (P(ii, 2) - z_low*h + corner(2))/(z_high*h + corner(2) - z_low*h + corner(2));
double w_z_high = (P(ii, 2) - z_low*h + corner(2))/(z_high*h + corner(2) - z_low*h + corner(2));

// x-weights
tripletList.push_back(T(ii, i000, w_x_low));
tripletList.push_back(T(ii, i010, w_x_low));
tripletList.push_back(T(ii, i001, w_x_low));
tripletList.push_back(T(ii, i011, w_x_low));

tripletList.push_back(T(ii, i100, w_x_high));
tripletList.push_back(T(ii, i110, w_x_high));
tripletList.push_back(T(ii, i101, w_x_high));
tripletList.push_back(T(ii, i111, w_x_high));

// y-weights
tripletList.push_back(T(ii + np, i000, w_y_low));
tripletList.push_back(T(ii + np, i001, w_y_low));
tripletList.push_back(T(ii + np, i101, w_y_low));
tripletList.push_back(T(ii + np, i100, w_y_low));

tripletList.push_back(T(ii + np, i010, w_y_high));
tripletList.push_back(T(ii + np, i110, w_y_high));
tripletList.push_back(T(ii + np, i011, w_y_high));
tripletList.push_back(T(ii + np, i111, w_y_high));

// z-weights
tripletList.push_back(T(ii + 2*np, i000, w_z_low));
tripletList.push_back(T(ii + 2*np, i100, w_z_low));
tripletList.push_back(T(ii + 2*np, i010, w_z_low));
tripletList.push_back(T(ii + 2*np, i110, w_z_low));

tripletList.push_back(T(ii + 2*np, i001, w_z_high));
tripletList.push_back(T(ii + 2*np, i101, w_z_high));
tripletList.push_back(T(ii + 2*np, i011, w_z_high));
tripletList.push_back(T(ii + 2*np, i111, w_z_high));


// ****** DEBUGGING ******
// std::cout << x_low << ", " << P(ii, 0) << ", " << x_high << std::endl;
// std::cout << y_low << ", " << P(ii, 1) << ", " << y_high << std::endl;
// std::cout << z_low << ", " << P(ii, 2) << ", " << z_high << std::endl;

// std::cout << (P(ii, 0) - x_low)/(x_high - x_low) << std::endl;
// std::cout << (P(ii, 1) - y_low)/(y_high - y_low) << std::endl;
// std::cout << (P(ii, 2) - z_low)/(z_high - z_low) << std::endl;
// ***********************

}

W.setFromTriplets(tripletList.begin(), tripletList.end());
W.finalize();

std::cout << "done." << std::endl;

////////////////////////////////////////////////////////////////////////////
}


134 changes: 133 additions & 1 deletion src/fd_partial_derivative.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include "fd_partial_derivative.h"
#include <iostream>

void fd_partial_derivative(
const int nx,
Expand All @@ -9,6 +10,137 @@ void fd_partial_derivative(
Eigen::SparseMatrix<double> & D)
{
////////////////////////////////////////////////////////////////////////////
// Add your code here

// Compute the derivatives for the desired direction
switch (dir)
{
case 1:
{
D.resize((nx-1)*ny*nz, nx*ny*nz);

// Using the triplet insertion method based on Eigen's documentation
typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;

for (int ii = 0; ii < nx*ny*nz; ii++)
{
int col = ii;
int col_other;

// Check if we're at an edge, in which case the finite difference is taken to be a forward difference
if (ii%nx == 0)
{
col_other = ii+1;
if (col >= (nx-1)*ny*nz)
col = (nx-1)*ny*nz - 1;
if (col_other >= nx*ny*nz)
col_other = nx*ny*nz - 1;

tripletList.push_back(T(col, col, -1.0));
tripletList.push_back(T(col, col_other, 1.0));
}
else
{
col_other = ii-1;
if (col_other >= (nx-1)*ny*nz)
col_other = (nx-1)*ny*nz - 1;
if (col_other < 0)
col_other = 0;

tripletList.push_back(T(col_other, col_other, -1.0));
tripletList.push_back(T(col_other, col, 1.0));
}
}

D.setFromTriplets(tripletList.begin(), tripletList.end());
D.finalize();
}
case 2:
{
D.resize(nx*(ny-1)*nz, nx*ny*nz);

// Using the triplet insertion method based on Eigen's documentation
typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;

for (int ii = 0; ii < nx*ny*nz; ii++)
{
int col = ii;
int col_other;

// Check if we're at an edge, in which case the finite difference is taken to be a forward difference
if (ii%(nx*ny) < nx)
{
col_other = ii+1*nx;
if (col >= nx*(ny-1)*nz)
col = nx*(ny-1)*nz - 1;
if (col_other >= nx*ny*nz)
col_other = nx*ny*nz - 1;

tripletList.push_back(T(col, col, -1.0));
tripletList.push_back(T(col, col_other, 1.0));
}
else
{
col_other = ii-1*nx;
if (col_other >= nx*(ny-1)*nz)
col_other = nx*(ny-1)*nz - 1;
if (col_other < 0)
col_other = 0;

tripletList.push_back(T(col_other, col_other, -1.0));
tripletList.push_back(T(col_other, col, 1.0));
}
}

D.setFromTriplets(tripletList.begin(), tripletList.end());
D.finalize();
}
case 3:
{
D.resize(nx*ny*(nz-1), nx*ny*nz);

// Using the triplet insertion method based on Eigen's documentation
typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;

for (int ii = 0; ii < nx*ny*nz; ii++)
{
int col = ii;
int col_other;

// Check if we're at an edge, in which case the finite difference is taken to be a forward difference
if (ii < nx)
{
col_other = ii+1*nx*ny;
if (col >= nx*ny*(nz-1))
col = nx*ny*(nz-1) - 1;
if (col_other >= nx*ny*nz)
col_other = nx*ny*nz - 1;

tripletList.push_back(T(col, col, -1.0));
tripletList.push_back(T(col, col_other, 1.0));
}
else
{
col_other = ii-1*nx*ny;
if (col_other >= nx*ny*(nz-1))
col_other = nx*ny*(nz-1) - 1;
if (col_other < 0)
col_other = 0;

tripletList.push_back(T(col_other, col_other, -1.0));
tripletList.push_back(T(col_other, col, 1.0));

}
}

D.setFromTriplets(tripletList.begin(), tripletList.end());
D.finalize();
}

};


////////////////////////////////////////////////////////////////////////////
}
Loading