Commit ba774134 authored by rtimothy's avatar rtimothy
Browse files

Merge branch 'dev' into 'master'

Dev

See merge request !1
parents 4eb6916c 73ee92bd
.vscode
test.py
tetgen/tetgen
TMSmesh/TMSmesh
TMSmesh/TMSmesh2
.DS_Store
*/.DS_Store
main/main
main/output.xml
*.xml
\ No newline at end of file
CXX= g++
CXXFLAGS= -std=c++11 -O3
main: main.cpp mem_gen.cpp box_gen.cpp
$(CXX) $(CXXFLAGS) -o ./main/$@ main.cpp
clean:
rm -f box_gen mem_gen ./main/main
.PHONY: all main clean
# Ion
# Ion Channel Mesh Generation Project
This is the public folder for the "Ion Channel Mesh Generation" team.
\ No newline at end of file
#### Required Packages
1. `fenics`
2. `tetgen`
3. `TMSmesh`
Make sure relative to project dir, tetgen's rel path is ./tetgen/tetgen. Similarly for TMSmesh. Make sure fenics is downloaded or active python3 install (venv or conda env or global install)
Links to each package's site below:
- `fenics`: https://fenicsproject.org/download/
- `tetgen`: https://www.wias-berlin.de/software/index.jsp?id=TetGen
- `TMSmesh`: http://lsec.cc.ac.cn/~lubz/Meshing.html
### Usage:
First, ensure that `fenics` has been properly installed by running the command `python3 -c "import fenics"`. If no errors are thrown, the installation was successful.
Next, ensure that the compiled `tetgen` binary has been named `tetgen` and moved into the folder `/tetgen/`.
Lastly, move the `TMSmesh` executable into the folder `/TMSmesh/` and set its executable permissions with `chmod +x <executable name>`. If on Linux, use the executable name `TMSmesh` and if on Mac, use the executable name `TMSmesh2`.
Makefile dependency tracking is faulty so make sure to `make clean` whenever updating source files
1. `make clean`
2. `make main`
3. Move desired protein .pqr file to `/pqr/` folder
4. `cd main`
5. Set configuration settings in config.txt
6. `./main`
7. Output files are `/main/membrane.pvd`, `/main/protein.pvd`, `/main/solvent.pvd`
### Configuration Settings
| Parameter | Type | Description |
| ------------- | ------------- | ------------- |
| input_file | <filename.(pqr or off)> | Input filename for the protein mesh. Can be in .pqr format (converted to .off with TMSmesh) or .off format |
| poly_file | <filename.poly> | Filename for intermediate .poly file of surface mesh, stored in `/poly/` and used by `tetgen` |
| surface_membrane_density | \<float> | Specifies the density of additional points inserted for the bottom and top of membrane mesh. |
| box_length(x) | \<int> | Specifies half box length in the x-dimension |
| box_width(y) | \<int> | Specifies half box width in y-dimension |
| top_box_height | \<int> | Specifies height of the box above membrane |
| bottom_box_height | \<int> | Specifies height of box below membrane |
| membrane_top_height | \<float> | Specifies z value for top of membrane |
| membrane_bottom_height | \<float> | Specifies z value for bottom of membrane |
| offset | \<int> | Specifies width of connective layer between membrane and top and bottom of the box mesh |
| tmsmesh_density_value | \<float> | Relate to the surface mesh density generated by TMSmesh. When this reduces,the density increases |
| tmsmesh_decay_rate | \<float> | Decay rate in the Gaussian surface used by TMSmesh |
| tmsmesh_isovalue | \<float> | Isovalue in the Gaussian surface and controls the volume enclosed by the Gaussian surface |
| tetgen_quality_q | \<float> | Specifies minimum radius-edge ratio in final volume mesh. Set to 0 for no quality constraint. |
| tetgen_quality_a | \<float> | Maximum volume constraint. No tetrahedra with volume larger than this value would be generated. Set to 0 for no volume constraint. |
| verbose | \<true or false> | Sets whether program is run verbose |
1 1 1.0000000000000000
3 1 -0.5000000000000000
5 1 0.3750000000000000
7 1 -0.3125000000000000
9 1 0.2734375000000000
11 1 -0.2460937500000000
13 1 0.2255859375000000
15 1 -0.2094726562500000
17 1 0.1963806152343750
19 1 -0.1854705810546875
21 1 0.1761970520019531
2 2 1.0000000000000000
4 2 -1.5000000000000000
6 2 1.8750000000000000
8 2 -2.1875000000000000
10 2 2.4609375000000000
12 2 -2.7070312500000000
14 2 2.9326171875000000
16 2 -3.1420898437500000
18 2 3.3384704589843750
20 2 -3.5239410400390625
3 3 1.5000000000000000
5 3 -3.7500000000000000
7 3 6.5625000000000000
9 3 -9.8437500000000000
11 3 13.5351562500000000
13 3 -17.5957031250000000
15 3 21.9946289062500000
17 3 -26.7077636718750000
19 3 31.7154693603515630
21 3 -37.0013809204101560
4 4 2.5000000000000000
6 4 -8.7500000000000000
8 4 19.6875000000000000
10 4 -36.0937500000000000
12 4 58.6523437500000000
14 4 -87.9785156250000000
16 4 124.6362304687500000
18 4 -169.1491699218750000
20 4 222.0082855224609400
5 5 4.3750000000000000
7 5 -19.6875000000000000
9 5 54.1406250000000000
11 5 -117.3046875000000000
13 5 219.9462890625000000
15 5 -373.9086914062500000
17 5 592.0220947265625000
19 5 -888.0331420898437500
21 5 1276.5476417541504000
6 6 7.8750000000000000
8 6 -43.3125000000000000
10 6 140.7656250000000000
12 6 -351.9140625000000000
14 6 747.8173828125000000
16 6 -1420.8530273437500000
18 6 2486.4927978515625000
20 6 -4084.9524536132812000
7 7 14.4375000000000000
9 7 -93.8437500000000000
11 7 351.9140625000000000
13 7 -997.0898437500000000
15 7 2368.0883789062500000
17 7 -4972.9855957031250000
19 7 9531.5557250976562000
21 7 -17020.6352233886720000
8 8 26.8125000000000000
10 8 -201.0937500000000000
12 8 854.6484375000000000
14 8 -2706.3867187500000000
16 8 7104.2651367187500000
18 8 -16339.8098144531250000
20 8 34041.2704467773440000
9 9 50.2734375000000000
11 9 -427.3242187500000000
13 9 2029.7900390625000000
15 9 -7104.2651367187500000
17 9 20424.7622680664060000
19 9 -51061.9056701660160000
21 9 114889.2877578735400000
10 10 94.9609375000000000
12 10 -902.1289062500000000
14 10 4736.1767578125000000
16 10 -18155.3442382812500000
18 10 56735.4507446289060000
20 10 -153185.7170104980500000
11 11 180.4257812500000000
13 11 -1894.4707031250000000
15 11 10893.2065429687500000
17 11 -45388.3605957031250000
19 11 153185.7170104980500000
21 11 -444238.5793304443400000
12 12 344.4492187500000000
14 12 -3961.1660156250000000
16 12 24757.2875976562500000
18 12 -111407.7941894531200000
20 12 403853.2539367675800000
13 13 660.1943359375000000
15 13 -8252.4291992187500000
17 13 55703.8970947265630000
19 13 -269235.5026245117200000
21 13 1043287.5726699828000000
14 14 1269.6044921875000000
16 14 -17139.6606445312500000
18 14 124262.5396728515600000
20 14 -642023.1216430664100000
15 15 2448.5229492187500000
17 15 -35503.5827636718750000
19 15 275152.7664184570300000
21 15 -1513340.2153015137000000
16 16 4733.8110351562500000
18 16 -73374.0710449218750000
20 16 605336.0861206054700000
17 17 9171.7588806152344000
19 17 -151334.0215301513700000
21 17 1324172.6883888245000000
18 18 17804.0025329589840000
20 18 -311570.0443267822300000
19 19 34618.8938140869140000
21 19 -640449.5355606079100000
20 20 67415.7405853271480000
21 21 131460.6941413879400000
#include <iostream>
#include <string>
#include <vector>
#include <map>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <limits>
#include <stdio.h>
#include "face.h"
using namespace std;
/**
* Get Protein Limits
*
* Search through protein off file to find minima and maxima for each axis
*
* @param protein_filename the filename of off protein file
* @param min_x ptr to int that willl contain min for x axis
* @param max_x ptr to in that will contain max for x axis
* @param min_y ptr to int that will contain min for y axis
* @param max_y ptr to int that will contain max for y axis
* @param min_z ptr to int that will contain min for z axis
* @param max_z ptr to int that will contain max for z axis
*
* return minima/maxima through ptrs passed in as arguments
**/
void getProteinLimits(string protein_filename, double *min_x, double *max_x, double *min_y, double *max_y, double *min_z, double *max_z) {
ifstream protein_file;
protein_file.open(protein_filename);
double minx = std::numeric_limits<double>::infinity();
double miny = std::numeric_limits<double>::infinity();
double minz = std::numeric_limits<double>::infinity();
double maxx = -std::numeric_limits<double>::infinity();
double maxy = -std::numeric_limits<double>::infinity();
double maxz = -std::numeric_limits<double>::infinity();
string trash;
protein_file >> trash; // Read "OFF" at the beginning of the file
int protein_nodes, protein_triangles, protein_edges;
protein_file >> protein_nodes >> protein_triangles >> protein_edges; // Read numbers of Nodes, Triangles and Edges of the file
string line;
getline(protein_file, line); // catches extra \n
// Iterate over only nodes and find min/max
for (int i = 0; i < protein_nodes; ++i) {
getline(protein_file, line);
istringstream pqrline(line);
double x, y, z;
pqrline >> x >> y >> z;
maxx = (x > maxx) ? x : maxx;
maxy = (y > maxy) ? y : maxy;
maxz = (z > maxz) ? z : maxz;
minx = (x < minx) ? x : minx;
miny = (y < miny) ? y : miny;
minz = (z < minz) ? z : minz;
}
protein_file.close();
// return min and max value
*min_x = minx;
*max_x = maxx;
*min_y = miny;
*max_y = maxy;
*min_z = minz;
*max_z = maxz;
}
/**
* MAIN FUNCTION OF box_gen.cpp (called in main.cpp)
*
* Generate the surface mesh (with the box) with output from TMSmesh
*
* @param config store map of config variables
* @param protein_filename the filename of protein surface mesh outputed by TMSmesh
**/
int boxGenMain(map<string,string>& config, string protein_filename) {
int option = 0; // 0 for .off, 1 for poly
const string output_file = "../poly/" + getVal(config, "poly_file");
if (output_file.find(".off") != string::npos) {
option = 0;
} else if (output_file.find(".poly") != string::npos) {
option = 1;
} else {
throw std::runtime_error("Invalid input file passed. Input file must be off or pqr.\n");
}
int n_dim = 3;
int n_nodes = 0, n_triangles = 0;
ostringstream node_out;
ostringstream triangle_out;
int delta = 1;
int nx = stoi(getVal(config, "box_length(x)")); // Complete length of box in x-axis is 2*nx
int ny = stoi(getVal(config, "box_width(y)")); // Complete length of box in y-axis is 2*ny
double mem_delta = 1 / stod(getVal(config, "surface_membrane_density")); // delta val for distance between additional nodes
int nz_bottom = stoi(getVal(config, "bottom_box_height")); // height of bottom z segment of box
int nz_top = stoi(getVal(config, "top_box_height")); // height of top z segment of box
int offset = stoi(getVal(config, "offset")); // height of connective layer
double membrane_bottom_height = 5; // protein_min_z
double membrane_top_height = 50; // protein_max_z
double protein_min_x = -66.1811;
double protein_max_x = 27.483;
double protein_min_y = -49.3957;
double protein_max_y = 49.1864;
getProteinLimits(protein_filename, &protein_min_x, &protein_max_x, &protein_min_y, &protein_max_y, &membrane_bottom_height, &membrane_top_height);
membrane_bottom_height = stod(getVal(config, "membrane_bottom_height"));
membrane_top_height = stod(getVal(config, "membrane_top_height"));
int nz_middle = int(membrane_top_height - membrane_bottom_height);
double shift_x = (-1) * nx;
double shift_y = (-1) * ny;
double shift_z = ((membrane_top_height + membrane_bottom_height) / 2) - ((nz_middle / 2) + (nz_bottom * delta * 2) + offset);
if (nz_top != nz_bottom) {
throw runtime_error("Different sizes for top and bottom layers of non-uniform faces not yet implemented.\n");
}
// Generate 6 faces one by one
Face xy_bottom(0, "z", 0, nx, ny, 0, 2 * delta);
xy_bottom.genNodes();
xy_bottom.shiftNodes(shift_x, shift_y, shift_z);
xy_bottom.genTriangles();
xy_bottom.genOutput(n_nodes, n_triangles, node_out, triangle_out, option);
int top_height = (nz_bottom * 2 * delta) + (nz_middle * delta) + (nz_top * 2 * delta) + (2 * offset);
Face xy_top(0, "z", top_height, nx, ny, 0, 2 * delta);
xy_top.genNodes();
xy_top.shiftNodes(shift_x, shift_y, shift_z);
xy_top.genTriangles();
xy_top.genOutput(n_nodes, n_triangles, node_out, triangle_out, option);
NonuniformFace yz_front("x", 0, nx, ny, nz_bottom, nz_middle, offset, delta);
yz_front.genNodes();
yz_front.shiftNodes(shift_x, shift_y, shift_z);
yz_front.genTriangles();
yz_front.genConnectorsBottom(ny);
yz_front.genConnectorsTop(ny);
yz_front.genOutput(n_nodes, n_triangles, node_out, triangle_out, option);
NonuniformFace yz_back("x", nx * 2 * delta, nx, ny, nz_bottom, nz_middle, offset, delta);
yz_back.genNodes();
yz_back.shiftNodes(shift_x, shift_y, shift_z);
yz_back.genTriangles();
yz_back.genConnectorsBottom(ny);
yz_back.genConnectorsTop(ny);
yz_back.genOutput(n_nodes, n_triangles, node_out, triangle_out, option);
NonuniformFace xz_left("y", 0, nx, ny, nz_bottom, nz_middle, offset, delta);
xz_left.genNodes();
xz_left.shiftNodes(shift_x, shift_y, shift_z);
xz_left.genTriangles();
xz_left.genConnectorsBottom(nx);
xz_left.genConnectorsTop(nx);
xz_left.genOutput(n_nodes, n_triangles, node_out, triangle_out, option);
NonuniformFace xz_right("y", ny * 2 * delta, nx, ny, nz_bottom, nz_middle, offset, delta);
xz_right.genNodes();
xz_right.shiftNodes(shift_x, shift_y, shift_z);
xz_right.genTriangles();
xz_right.genConnectorsBottom(nx);
xz_right.genConnectorsTop(nx);
xz_right.genOutput(n_nodes, n_triangles, node_out, triangle_out, option);
// File output
ofstream poly_file;
poly_file.open(output_file);
ifstream protein_file;
protein_file.open(protein_filename);
string trash;
protein_file >> trash; // Read "OFF" at the beginning of the file
int protein_nodes, protein_triangles, protein_edges;
protein_file >> protein_nodes >> protein_triangles >> protein_edges;
string line;
getline(protein_file, line); // Catches extra \n
if (!option) {
// .off output
for (int i = 0; i < protein_nodes; ++i) {
getline(protein_file, line);
node_out << line << "\n";
}
for (int i = 0; i < protein_triangles; ++i) {
string line;
getline(protein_file, line);
istringstream triangle(line);
int n1, n2, n3, trash;
triangle >> trash >> n1 >> n2 >> n3;
n1 += n_nodes;
n2 += n_nodes;
n3 += n_nodes;
triangle_out << "3\t" << n1 << "\t" << n2 << "\t" << n3 << "\n";
}
protein_file.close();
cout << "OFF\n";
cout << n_nodes + protein_nodes << "\t" << n_triangles + protein_triangles << "\t0\n";
cout << node_out.str();
cout << triangle_out.str();
} else {
// .poly output
for (int i = 0; i < protein_nodes; ++i) {
getline(protein_file, line);
node_out << (n_nodes + i + 1) << "\t" << line << "\n";
}
for (int i = 0; i < protein_triangles; ++i) {
string line;
getline(protein_file, line);
istringstream triangle(line);
int n1, n2, n3, trash;
triangle >> trash >> n1 >> n2 >> n3;
n1 += n_nodes;
n2 += n_nodes;
n3 += n_nodes;
triangle_out << "1\n";
triangle_out << "3\t" << (n1 + 1) << "\t" << (n2 + 1) << "\t" << (n3 + 1) << "\n";
}
// Create additional nodes
vector<Node> add_nodes;
for (double i = -nx + mem_delta; i <= (nx - mem_delta); i += mem_delta) {
for (double j = -ny + mem_delta; j <= (ny - mem_delta); j += mem_delta) {
if (!((i >= protein_min_x && i <= protein_max_x) &&
(j >= protein_min_y && j <= protein_max_y))) {
add_nodes.emplace_back(i, j, membrane_bottom_height);
}
}
}
for (double i = -nx + mem_delta; i <= (nx - mem_delta); i += mem_delta) {
for (double j = -ny + mem_delta; j <= (ny - mem_delta); j += mem_delta) {
if (!((i >= protein_min_x && i <= protein_max_x) &&
(j >= protein_min_y && j <= protein_max_y))) {
add_nodes.emplace_back(i, j, membrane_top_height);
}
}
}
for (int i = 0; i < add_nodes.size(); ++i) {
node_out << (n_nodes + i + 1) << "\t" << add_nodes[i].x << "\t" << add_nodes[i].y << "\t" << add_nodes[i].z << "\n";
}
n_nodes += add_nodes.size();
// File output
poly_file << "#POLY FILE OUTPUT\n";
poly_file << (n_nodes + protein_nodes) << "\t3\t0\t0\n";
poly_file << node_out.str();
poly_file << (6 + protein_triangles) << "\t0\n";
poly_file << triangle_out.str();
poly_file << "# number of holes\n";
poly_file << "0\n";
poly_file << "# number of regions\n";
poly_file << "0\n";
}
poly_file.close();
return 0;
}
#include <iostream>
#include <vector>
#include <string>
#include <sstream>
using namespace std;
extern string getVal(map<string,string>& config, string key);
struct Triangle {
int node1, node2, node3;
Triangle(int node1_in, int node2_in, int node3_in) : node1(node1_in), node2(node2_in), node3(node3_in) {}
Triangle() : node1(0), node2(0), node3(0) {}
};
struct Node {
double x, y, z;
Node(double x_in, double y_in , double z_in) : x(x_in), y(y_in), z(z_in) {}
Node() : x(0.0), y(0.0), z(0.0) {}
};
class Face {
public:
Face(int offset_in, string fixed_axis_in, int fixed_axis_val_in, int x_dim_in, int y_dim_in, int z_dim_in, int delta_in) :
offset(offset_in), fixed_axis(fixed_axis_in), fixed_axis_val(fixed_axis_val_in),
x_dim_num(x_dim_in), y_dim_num(y_dim_in), z_dim_num(z_dim_in), delta(delta_in),
x_dim(x_dim_num* delta), y_dim(y_dim_num* delta), z_dim(z_dim_num* delta)
{}
inline size_t numNodes() { return nodes.size(); }
inline size_t numTriangles() { return triangles.size(); }
/**
* Generate all nodes for a face in a uniform grid
**/
void genNodes() {
int i_lim, j_lim;
setLoopLimit(i_lim, j_lim, true);
bool xIsFixed = fixed_axis.compare("x") == 0;
bool yIsFixed = fixed_axis.compare("y") == 0;
bool zIsFixed = fixed_axis.compare("z") == 0;
for (int i = 0; i <= i_lim; i += delta) {
for (int j = 0; j <= j_lim; j += delta) {
if (xIsFixed) {
nodes.emplace_back(fixed_axis_val, i, j + offset);
}
else if (yIsFixed) {
nodes.emplace_back(i, fixed_axis_val, j + offset);
}
else {
nodes.emplace_back(i, j + offset, fixed_axis_val);
}
}
}
}
/**
* Shift each node in face by deltaX, deltaY and deltaZ
* @param delta_x shift in X dimension
* @param delta_y shift in Y dimension
* @param delta_z shift in Z dimension
**/
void shiftNodes(double delta_x, double delta_y, double delta_z) {
for (Node & n : nodes) {
n.x += delta_x;
n.y += delta_y;
n.z += delta_z;
}
}
/**
* Generate triangles by connecting nodes in a face
**/
void genTriangles() {
int i_lim, j_lim;
setLoopLimit(i_lim, j_lim, false);
for (int i = 0; i < i_lim; ++i) {
for (int j = 0; j < j_lim; ++j) {
triangles.emplace_back(i * (j_lim + 1) + j, (i * (j_lim + 1)) + (j + 1), (i + 1) * (j_lim + 1) + j);
triangles.emplace_back((i * (j_lim + 1)) + (j + 1), (i + 1) * (j_lim + 1) + j, (i + 1) * (j_lim + 1) + (j + 1));
}
}
}
/**
* Helper function for genOutput
* Prints triangles and nodes to respective output streams
*
* @param n_nodes tracks number of nodes output
* @param n_triangles tracks number of triangles output
* @param node_out node output stream
* @param triangle_out triangle output stream
* @param option specifies file format (0 for .off, 1 or 2 for .poly)
* 1 for a uniform face, 2 if part of a nonuniform face
**/
void genOutput(int& n_nodes, int& n_triangles, ostringstream& node_out, ostringstream& triangle_out, int option) {
if(option == 0){
output_helper(n_nodes, n_triangles, node_out, triangle_out, 0);
}
else if(option == 1){
triangle_out << "# a new facet begins here \n";
triangle_out << triangles.size() << "\n";
output_helper(n_nodes, n_triangles, node_out, triangle_out, 1);
}
else{
// triangle_out << "# a facet continues here \n";
output_helper(n_nodes, n_triangles, node_out, triangle_out, 1);
}
}
private:
vector<Node> nodes;
vector<Triangle> triangles;
string fixed_axis;
int fixed_axis_val;
int offset;
int delta;
int x_dim_num, y_dim_num, z_dim_num;
int x_dim, y_dim, z_dim;
/**
* Helper function for genOutput
* Prints triangles and nodes to respective output streams
*
* @param n_nodes tracks number of nodes output
* @param n_triangles tracks number of triangles output
* @param node_out node output stream
* @param triangle_out triangle output stream
* @param offset specifies offset for node and triangle indices when output
**/
void output_helper(int& n_nodes, int& n_triangles, ostringstream& node_out, ostringstream& triangle_out, int offset){
for (int i = 0; i < triangles.size(); ++i) {
triangle_out << "3\t" << triangles[i].node1 + n_nodes + offset
<< "\t" << triangles[i].node2 + n_nodes + offset
<< "\t" << triangles[i].node3 + n_nodes + offset << "\n"